2. Generalized Pólya point process with stationary increments: $\Psi$-GPP

The Pólya point process has been extended to the generalized Pólya process (GPP) [Reference Cha7, Reference Konno27]. A GPP $N \equiv \{N(t)\,{:}\, t \ge 0\}$ is a Markov point process with stochastic intensity (defined in terms of the internal histories $\mathcal{H}_{t}$; see, e.g., [Reference Baccelli and Bremaud3, Section 1.8]) by

(1) \begin{equation}\lambda (t) \equiv \lambda (t \mid \mathcal{H}_{t}) \equiv (\gamma N(t-) + \beta)\kappa (t),\end{equation}

where $N(0) = 0$, $\gamma$ and $\beta$ are positive constants, $\kappa (t)$ is a positive integrable real-valued function and $\equiv$ denotes equality by definition. The classical Pólya point process is the special case of (1) with $\beta = 1$ and

(2) \begin{equation}\kappa (t) = \frac{1}{\gamma t + 1}, \qquad t \ge 0.\end{equation}

Many properties of the GPP were deduced in [Reference Cha7] by exploiting the restarting property.

Theorem 1 of [Reference Cha7] establishes the joint distribution of a GPP by exploiting the restarting property. As a consequence, the marginal distribution of a GPP starting at $N(0) = 0$ has a simple form.

Proposition 2. (Negative binomial marginal distribution [Reference Cha7]) If N is a GPP with parameter triple $(\kappa (t), \gamma, \beta)$, then N(t) has a negative binomial distribution with probability mass function $\mathbb{P}(N(t) = k) \equiv f(k;\ r, p(t)) = C(\beta, \gamma, k) (1-p(t))^r p(t)^k$, $k = 0,1, 2, \ldots$, where $r = \beta/\gamma$, $p(t) = 1 - \exp\!\{-\gamma K (t)\}$, $K(t) \equiv \int_{0}^{t} \kappa (s) \, \mathrm{d} s$, $t \ge 0$, and $C(\beta, \gamma, k) = \Gamma ((\beta/\gamma) + k)/\Gamma ((\beta/\gamma)) k!$, with $\Gamma$ being the gamma function, so N(t) has mean and variance

\begin{displaymath}\mathbb{E}[N(t)] = \frac{r p(t)}{1-p(t)} , \quad \mathrm{Var}(N(t)) = \frac{r p(t)}{(1-p(t))^2}, \qquad t \ge 0.\end{displaymath}

For a general function $\kappa (t)$ the time-varying behavior can be complicated, but it simplifies for the classical Pólya point process and closely related processes (allowing $\beta \not= 1$). Indeed, for $\kappa (t)$ in (2) the GPP is a (strictly) stationary point process, i.e. the joint distribution of any k increments is independent of time shifts, as we show next. In the spirit of [Reference Fendick16], we thus call the GPP with triple $(\kappa (t), \gamma, \beta)$ for $\kappa (t)$ in (2) a $(\beta,\gamma)$ $\Psi$-GPP (again, with $\Psi$ as a mnemonic for SI).

Theorem 1. (A stationary point process: $\Psi$-GPP) Consider a GPP with parameter triple $(\kappa (t), \gamma, \beta)$. If $\kappa (t)$ is given by (2), then $1 - p(t) = \kappa (t)$,

(3) \begin{equation}\mathbb{E}[N(t)] = \beta t,\quad \mathrm{Var}(N(t)) = \beta t(1 + \gamma t), \qquad t \ge 0.\end{equation}

Moreover, the joint distribution of k increments, $N(s_i+t_i+h) - N(s_i+h)$, $1 \le i \le k$, is independent of $h > 0$ for all h, so that N is a stationary stochastic point process with

(4) \begin{equation}\mathrm{Cov}(N(s), N(t)) = \beta s(1 + \gamma t), \qquad 0 \le s \le t < \infty.\end{equation}

Proof. For $\kappa$ in (2),

\begin{displaymath}K(t) \equiv \int_{0}^{t} (\gamma s + 1)^{-1} \, \mathrm{d} s = \gamma^{-1}\log{(\gamma t + 1)}, \qquad t \ge 0,\end{displaymath}

so that $1 - p(t) = \mathrm{e}^{-\gamma K(t)} = \kappa (t)$. Then the conclusion about the distribution of a single increment follows from the displayed distribution of an increment in [Reference Cha7, Theorem 1(ii)]. For the covariance in (4), write $\mathrm{Var}(N(t-s)) = \mathrm{Var}(N(t)-N(s)) = \mathrm{Var}(N(t)) + \mathrm{Var}(N(s))-2 \mathrm{Cov}(N(s),N(t))$, where the first equality follows from stationary increments, to obtain $\mathrm{Cov}(N(s), N(t)) = [\mathrm{Var}(N(t)) + \mathrm{Var}(N(s)) - \mathrm{Var}(N(t-s))]/2$, and then use the variance formula in (3).

For the joint distribution of k increments, we first observe that, without loss of generality, we can assume that the k increments are disjoint and contiguous, so that they represent a partition of a fixed interval $(s,s+t]$ into finitely many subintervals. (We are thus initially adding intervals and determining the joint distribution for more intervals.) We then apply [Reference Cha7, Theorem 3 and Remark 3] to conclude that the conditional distribution for the sequence of times when N increases on $(s,s+t]$ given that $N(s + t) - N(s)=k$ is, first, independent of s, and second, is itself the same as that of the order statistics of k i.i.d. random variables, each with probability density function

\begin{equation*}f(x)\equiv \frac{\gamma \kappa (x) \exp\!\left(\gamma K(x)\right)}{\exp\!\left(\gamma K(t)\right)-1}, \qquad 0\le x\le t,\end{equation*}

assuming $s = 0$. If $\kappa (t)$ is given by (2), then $f(x)=1/t$; i.e. the probability density function (PDF) of the uniform distribution. That in turn implies that the joint distribution of the k disjoint and contiguous increments $N(s_i+t_i+h) - N(s_i+h)$, $1 \le i \le k$, all within some larger interval $(s+h, s + t + h)$, is independent of $h > 0$ for all h. That follows because the conditioning event that $N(s + t+ h) - N(s + h) =k$ has a distribution that is independent of h, and then the conditional distribution of the points within the interval given that $N(s + t+ h) - N(s + h) =k$ is also independent of h. Finally, we go from the case of contiguous intervals to the original case by integrating out in order to obtain the desired marginal distribution; e.g. when we let the probability that one interval is $\le \infty$, then that interval drops out.

We now show that any GPP with general $\kappa (t)$ in (1) can be expressed as a unique deterministic time-transformation of a $\Psi$-GPP based on (2). As a consequence, limits established for $\Psi$-GPPs will extend to GPPs by the continuous mapping theorem. Let $\stackrel{\rm d}{=}$ denote equality in distribution, including for stochastic processes.

Theorem 2. (Characterization of GPP) Let $\hat{N}$ be the GPP with parameter triple $(\kappa (t), \gamma, \beta)$ satisfying (1); let N be the $\Psi$-GPP with parameter triple $(\kappa (t), \gamma, \beta)$ satisfying (1) and (2); and let

(5) \begin{equation}M(t) \equiv \int_{0}^{t} \mu (s) \, \mathrm{d} s, \qquad t \ge 0,\end{equation}

where $\mu$ is a continuous positive function on $[0, \infty)$. Then $\{\hat{N}(t)\,{:}\, t \ge 0\} \stackrel{\rm d}{=} \{N(M(t))\,{:}\, t \ge 0\}$ if and only if

(6) \begin{equation}M(t) = \gamma^{-1} (\mathrm{e}^{\gamma K(t)} -1), \quad t \ge 0, \qquad \mbox{where} \ K(t) \equiv \int_{0}^{t} \kappa (s) \, \mathrm{d} s.\end{equation}

Proof. The stochastic intensity of N can be interpreted as

\begin{displaymath}\lambda (t) \equiv \lambda (t \mid \mathcal{H}_{t-}) \equiv \lim_{h \rightarrow 0}{\frac{\mathbb{P}(N(t+h) - N(t) = 1 \mid \mathcal{H}_{t-})}{h}},\end{displaymath}

where $\{\mathcal{H}_t\,{:}\, t \ge 0\}$ is the family of histories associated with N. Let $C(t) = N(M(t))$, $t \ge 0$, be the composition process and let $\lambda_{C} (t)$ be its stochastic intensity function. Let $\{\mathcal{C}_t\,{:}\, t > 0\}$ be the histories associated with C. Because of (5), $\mathcal{C}_t = \mathcal{N}_{M(t)}$, $t \ge 0$. Then,

\begin{align*}\lambda_{C} (t) & \equiv \lambda_{C} (t \mid \mathcal{C}_{t-})\equiv \lim_{h \rightarrow 0}{\frac{\mathbb{P}(C(t+h) - C(t) = 1 \mid \mathcal{C}_{t-})}{h}} \nonumber \\[3pt]& = \lim_{h \rightarrow 0}{\frac{\mathbb{P}(N(M(t+h)) - N(M(t)) = 1 \mid \mathcal{N}_{M(t)-})}{h}} \nonumber \\[3pt]& = \lim_{h \rightarrow 0}{\frac{\mathbb{P}(N(M(t) + \delta(h))) - N(M(t)) = 1 \mid \mathcal{N}_{M(t)-}))}{\delta(h)}\frac{\delta(h)}{h}},\end{align*}

where $\delta \equiv \delta(h) \equiv M(t + h) - M(t) = \mu (t) h + o(h)$ as $h \downarrow 0$. Hence,

\begin{equation*}\lambda_{C} (t) = \lambda (M(t)) \mu (t) = \frac{(\gamma N((M(t)) + \beta) \mu (t)}{\gamma M(t) + 1} ,\end{equation*}

where we use the specific form of $\kappa (t)$ in (2) because N is a $\Psi$-GPP.

Now let $\hat{\lambda} (t)$ be the stochastic intensity associated with $\hat{N}$. Then $\hat{\lambda} (t) = (\gamma \hat{N} (t) + \beta) \kappa (t)$, $t \ge 0$. Under the condition that $\hat{N} (t) = N(M(t))$, we have $\gamma \hat{N} (t) + \beta = \gamma N((M(t)) + \beta$. Hence, $\lambda_{C} (t) = \hat{\lambda} (t)$ for all $t \ge 0$ if and only if

(7) \begin{equation}\frac{\mu (t)}{\gamma M(t) + 1} = \kappa (t), \qquad t \ge 0.\end{equation}

Integrating both sides of (7), we get (6), as claimed. In closing, we mention two alternative proofs. The first can be based on the martingale characterization of the stochastic intensity as in [Reference Baccelli and Bremaud3, Section 1.8]. Then $\lambda (t) \equiv \lambda (t \mid \mathcal{H}_t)$ is the stochastic intensity of N in (1) if

\begin{displaymath}\mathbb{E}[N(b) - N(a) \mid \mathcal{H}_a] = \mathbb{E}\bigg[\int_{a}^{b} \lambda (s) \, \mathrm{d} s \mid \mathcal{H}_a\bigg]\end{displaymath}

for all intervals (a, b]. The second is to directly show that the finite-dimensional distributions are the same.

Remark 1. (Index of dispersion) To better understand the impact of the variability as a function of time in an arrival process on the performance of a queueing model, we have shown in [Reference Fendick and Whitt17, Reference Whitt and You.38, Reference Whitt and You.39] that it is often helpful to look at the index of dispersion for counts, which for a $(\beta,\gamma)$ $\Psi$-GPP is

\begin{displaymath} I(t) \equiv \frac{\mathrm{Var}(N(t))}{EN(t)} = 1 + \gamma t, \qquad t \ge 0.\end{displaymath}

From this equation, we see that the variability increases without bound as t increases by this measure, consistent with Theorem 3. The index of dispersion for counts is also considered in [Reference Mirtchev31, Reference Mirtchev and Goleva32, Reference Mirtchev and Ganchev33] under the name ‘peakedness’, which is often used to describe traffic variability, but more commonly in a different way; see [Reference Li and Whitt28] and references therein.

We next apply Theorem 1 to characterize the nature of the path-dependent behavior for a $(\beta,\gamma)$ $\Psi$-GPP. We do that by the following non-ergodic LLN; see, e.g., [Reference Gregoire20, Section 5.1] and references therein.

Proof. Because the increments $N(n) - N(n-1)$, $n \ge 1$, form a stationary sequence, we can apply the Birkhoff ergodic theorem as in [Reference Breiman6, Theorem 6.2.1] to establish the almost sure convergence of $N(n)/n$ as $n \rightarrow \infty$. The almost sure convergence of $N(t)/t$ as $t \rightarrow \infty$ is an easy consequence, e.g. using $N(\lfloor nt \rfloor) \le N(t) \le N(\lfloor nt \rfloor + 1)$ for $\lfloor nt \rfloor \le t < \lfloor nt \rfloor + 1$. Next, from (3), we see that $\mathbb{E}[N(t)/t] = \beta$ for all $t > 0$ and $\mathrm{Var}(N(t)/t) = \beta(1 + \gamma t)/t \rightarrow \beta \gamma$ as $t \rightarrow \infty$. The limiting gamma distribution for $N(t)/t$ is obtained by taking the limit as $t \rightarrow \infty$ of the characteristic function $\phi_{N(t)/t} (s) = \phi_{N(t)} (s/t)$. In particular, using Taylor-series asymptotics in the last step, we obtain

\begin{eqnarray*} \phi_{t^{-1} N(t)} (s) & = & \left(\frac{1 - p(t)}{1 - p(t)\mathrm{e}^{is/t}}\right)^{\beta/\gamma} = \left(\frac{1}{1 + [p(t)/(1 - p(t))](1 - \mathrm{e}^{is/t})}\right)^{\beta/\gamma} \nonumber \\ & = & \left(\frac{1}{1 + \gamma t (1 - \mathrm{e}^{is/t})}\right)^{\beta/\gamma} = \left(\frac{1}{1 + \gamma t (1 - (1 + (is/t) + O(1/t^2))}\right)^{\beta/\gamma} \nonumber \\ & \rightarrow & \left(\frac{1}{1 - is \gamma}\right)^{\beta/\gamma} \quad\mbox{as}\quad t \rightarrow \infty.\end{eqnarray*}

Finally, we recognize the limit as the characteristic function of the claimed gamma distribution. Convergence of characteristic functions then applies the convergence in distribution by [Reference Feller15, Theorem XV.2].

Theorem 3 implies that the pure birth process N has a limiting rate as $t \rightarrow \infty$, but that rate is random.

Proof. Multiply and divide by t in (1) and observe that the numerator converges to $\gamma L$ by Theorem 3, while the denominator converges to $\gamma$.

We illustrate Theorem 3 by showing the results of a simulation experiment. Figure 1 plots the empirical distribution of $N(100)/100$ based on 50,000 i.i.d. samples (left) and 25 individual sample paths of $N(t)/t$ over [0, 200] for a stationary Pólya point process with parameter pair $(\beta, \gamma) = (1, 1)$, assuming that (2) holds (right).

Figure 1: Display of the empirical distribution of $N(100)/100$ (left) and 25 individual sample paths of $N(t)/t$ over [0, 200] (right) for a stationary Pólya point process N with $(\beta, \gamma) = (1, 1)$.

Figure 1 shows that much of the uncertainty about the long-run average disappears after an initial time t, e.g. $t = 50$, but the process $N_u (t) \equiv N(t + u) - N(u)$, $t \ge 0$, given that $N(u) = m$, is itself a new GPP by the restart property in Proposition 1. Additional results about this conditional process are given in [Reference Cha7]. We now apply [Reference Cha7] to further quantify the impact of starting at time u. The following result is elementary to prove, but the relations are revealing.

Corollary 2. (Impact of starting later) If $N_u (t) \equiv N(u + t) - N(u)$ conditioned on $N(u) = m$, where N is the $\Psi$-GPP with parameter triple $(\kappa (t),\beta, \gamma)$ for $\kappa (t)$ in (2), then $N_u$ is a new $\Psi$-GPP with parameter triple $(\kappa_{u} (t),\beta_{u}, \gamma_{u})$ for $\kappa_{u} (t) = \kappa (u + t)$, $\beta_u \equiv \beta + m \gamma$ and $\gamma_u \equiv \gamma$. Consequently, $N_{u} (t)$ has a negative binomial distribution for each $t \ge 0$ with

\begin{displaymath}\mathbb{E}[N_u (t)] = \left(\frac{\beta_u}{\gamma_{u}}\right)\left(\frac{p_{u} (t)}{1 - p_{u} (t)}\right)= \frac{(\beta + m \gamma)t}{\gamma u + 1},\end{displaymath}

where

\begin{displaymath}p_{u} (t) \equiv \frac{\eta_{u} t}{\eta_{u} t + 1}, \qquad \eta_{u} \equiv \frac{\gamma}{\gamma u + 1},\end{displaymath}

so that $\mathbb{E}[N_u (t)] \le \mathbb{E}[N (t)]$ if and only if $(m/u) \le \gamma$. Moreover, for $0 \le s \le t$,

\begin{align*} \mathrm{Cov}(N_{u} (s), N_{u} (t)) & = (\beta_{u}\eta_u/\gamma_u) s(1 + \eta_{u} t) = \mathbb{E}[N_{u} (s)](1 + \eta_{u} t) , \nonumber \\[3pt] I_{u} (t) & \equiv \frac{\mathrm{Var}(N_u (t))}{\mathbb{E}[N_{u} (t)]} = 1 + \eta_{u} t,\end{align*}

where $\eta_{u}$ is strictly decreasing in u with $\eta_{0} = \gamma$ and $\eta_{u} \rightarrow 1$ as $u \rightarrow \infty$.

Proof. This result mostly follows from the restart property from [Reference Cha7] stated here in Proposition 1. We also apply the explicit distribution of $N_{u} (t)$ given in [Reference Cha7, (20)]. The stationarity property follows from essentially the same proof as for Theorem 1.

We close this section by comparing the $\Psi$-GPP to another self-exciting process.

Remark 3. (Comparison to the Hawkes process.) The nondegenerate limit in Theorem 2 makes the $\Psi$-GPP quite different from the widely applied Hawkes [Reference Hawkes23, Reference Hawkes24] process and most of its variants. For the basic Hawkes process in [Reference Hawkes23, (8)], instead of (1) we have

\begin{displaymath} \lambda (t) \equiv \lambda (t \mid \mathcal{H}_{t-}) \equiv \nu + \int_{-\infty}^{t-} g(t-u) \, \mathrm{d} N(u) , \end{displaymath}

where N is defined over the entire real line and g is a nonnegative kernel satisfying $\eta \equiv \int_{0}^{\infty} g(u) \, \mathrm{d} u < 1$, so that the stationary rate is $\nu/(1 - \eta)$. Hawkes processes are alternative self-exciting processes, but they are stationary and ergodic point processes [Reference Bacray, Delattre, Hoffman and Muzy4]. Thus, Hawkes processes are not path dependent. For applications of Hawkes processes to queues, see [Reference Daw and Pender9, Reference Gao and Zhu18].

3. Convergence to a $\Psi$-GMP

We now show that a properly scaled sequence of the superpositon of i.i.d. $\Psi$-GPPs in Section 2 converges to a $\Psi$-GMP, the Gaussian Markov process with stationary increments studied in [Reference Fendick16]. (In fact, [Reference Fendick16] focuses on a multivariate $\Psi$-GMP.) We obtain all possible univariate $\Psi$-GMPs exhibiting positive dependence, as we explain in Section 8.

For $n \ge 1$, let

(8) \begin{equation}A^n = N^{1} + \cdots + N^{n}\end{equation}

be the sum of n i.i.d. GPPs each with parameter triple $(\kappa (t), \gamma, \beta)$. We first note that our superposition process is another GPP.

We now apply the usual FCLT spatial scaling, but without scaling time by n (as in [Reference Whitt36, (2.1) on p. 226 or (8.4) on p. 320]). In particular, for $n \ge 1$, let

(9) \begin{equation}A_n (t) \equiv n^{-1/2}(A^n (t) - \beta nt), \qquad t \ge 0 . \end{equation}

Let $\Rightarrow$ denote convergence in distribution and let $D \equiv D[0, \infty)$ be the usual function space of right-continuous real-valued functions [Reference Billingsley5, Reference Whitt36].

Theorem 4. (FCLT for the superposition process.) Consider the scaled superposition process $A_n (t)$ in (9). For $\kappa (t)$ in (2), so that in (8) $N^1$ is a $(\beta, \gamma)$ $\Psi$-GPP while $A^n$ is an $(n\beta, \gamma)$ $\Psi$-GPP,

(10) \begin{equation}A_n \Rightarrow A \text{ in } D \qquad \text{as } n \rightarrow \infty,\end{equation}

where A is a $\Psi$-GMP, i.e. a zero-mean Gaussian Markov process with stationary increments and covariance function

(11) \begin{equation}\mathrm{Cov}(A(s), A(t)) = \mathbb{E}[A(s)A(t)] = \beta s (1 + \gamma t) = \mathrm{Cov}(N^1 (s), N^1 (t)).\end{equation}

The limit A also satisfies the stochastic differential equation

(12) \begin{equation}\mathrm{d} A(t) = \mu (t) A(t) \, \mathrm{d} t + \sigma \, \mathrm{d} B(t), \qquad t \ge 0,\end{equation}

where $A(0) \equiv 0$, B is standard Brownian motion,

(13) \begin{equation}\mu(t) \equiv \frac{\beta\gamma}{\beta + \beta \gamma t} = (t + (1/\gamma))^{-1} , \qquad \sigma = \sqrt{\beta}.\end{equation}

Proof. For the limit in (10), we apply Hahn’s [Reference Hahn21] FCLT for sums of processes in [Reference Whitt36, Theorem 7.2.1]. We verify the moment inequality conditions in that theorem in Section 7. The stochastic differential equation characterization in (12) and (13) follows from [Reference Fendick16, Theorem 3]. A Gaussian process with that covariance kernel is a Markov process by [Reference Doob12, Theorem 8.1, p. 233].

4. Convergence to a $\Psi$-GMP with drift

For stable queueing models, there tends to be a negative drift in the potential net input process. Hence, in this section we consider a modification of the FCLT in Theorem 4 to produce a drift in the $\Psi$-GMP limit process. For that purpose, let $\{\mu_n\,{:}\, n \ge 1\}$ be a sequence of real numbers that satisfies

(14) \begin{equation}\mu_n \rightarrow 1 \quad\mbox{and}\quad \sqrt{n}(\mu_n - 1) \rightarrow \mu \qquad \text{as } n \rightarrow \infty.\end{equation}

We are primarily interested in the case $\mu < 0$. Let $A^{d,n} (t) = A^n (\mu_n t)$ and

(15) \begin{equation}A^{d}_n (t) = n^{-1/2}(A^{d,n} (t) - \beta nt) = n^{-1/2}(A^n (\mu_n t) - \beta nt), \qquad t \ge0.\end{equation}

Let $e \equiv e(t) = t$, $t \ge 0$, be the identity function in D, and let $D^k$ be the usual k-fold product space.

Proof. We apply the continuous mapping argument for composition with centering as in [Reference Whitt36, Section 13.3]. For that purpose, let $(\bar{M}_n (t), M_n (t)) \equiv \left(\mu_n t, \sqrt{n}(\mu_n - 1)t\right)$. Note that $A^{d}_n = A_n \circ \bar{M}_n + \beta M_n$, where $\circ$ denotes the composition map. It is elementary that $(\bar{M}_n , M_n ) \rightarrow (e, \mu e)$ in $D^2$ as $n \rightarrow \infty$. Then apply [Reference Whitt36, Theorem 11.4.5] with Theorem 4 above to get the joint convergence $(A_n, \bar{M}_n , M_n) \Rightarrow (A, e, \mu e)$ in $D^3$. Then the limit preservation in [Reference Whitt36, Theorem 13.3.1] yields $A^{d}_n = (A_n \circ \bar{M}_n + \beta M_n) \Rightarrow A + \mu \beta e$ in D as $n \rightarrow \infty$.

We now state three properties of a GMP with drift. The first two provide additional characterization of the path-dependent behavior.

Proposition 4 shows that conditioning on the history induces the process to have a new constant drift.

Let $\mathrm{Cor}(X,Y) \equiv \mathrm{Cov}(X,Y)/\sqrt{\mathrm{Var}(X)\mathrm{Var}(Y)}$ be the correlation function.

Corollary 4. (Correlation between non-overlapping time intervals, Proposition 2 of [Reference Fendick16].) For $A^d$ as in Proposition‘ 4,

\begin{displaymath}\mathrm{Cor}\big(A^d (t + s + u) - A^d (t + u), A^d (t + s) - A^d (t)\big) = \frac{\gamma s}{\gamma s + 1}\end{displaymath}

for all $t \ge 0$ and $u \ge s \ge 0$.

Corollary 4 concludes that the correlation between increments over non-overlapping intervals of equal lengths depends on the length of the intervals (s here) but not at all on the separation between the intervals (u here). The following corollary gives the limiting distribution. Let $\Phi (x) \equiv \mathbb{P}(N(0,1) \le x)$ be the standard normal cumulative distribution function (CDF) and let $\Phi^c (x) \equiv 1 - \Phi (x)$.

Corollary 5. (Asymptotics for the CDF.) For $A^d$ as in Proposition 4,

\begin{displaymath} \lim_{t \rightarrow \infty}{\mathbb{P}(A^d (t) \le x)} = \Phi(-\omega/\sqrt{\beta\gamma})\end{displaymath}

for all x and all $t \ge 0$. Moreover,

(16) \begin{align} \begin{split} \lim_{x \rightarrow -\infty} \lim_{t \rightarrow \infty}{\mathbb{P}(A^d (t) \le x)} & = \Phi(-\omega/\sqrt{\beta\gamma}) , \\ \lim_{x \rightarrow +\infty} \lim_{t \rightarrow \infty}{\mathbb{P}(A^d (t) > x)} & = \Phi^{c}(-\omega/\sqrt{\beta\gamma}). \end{split}\end{align}

Proof. To get the first result we can directly take the limit in the Gaussian distribution of $A^d (t)$, which is

\begin{displaymath} \mathbb{P}(A^d (t) \le x) = \Phi\left(\frac{(x - \omega t)}{\sqrt{t(\beta + \beta\gamma t)}}\right)\!, \end{displaymath}

from which (16) follows directly.

Corollary 5 can be understood by recognizing that the standard deviation is the same order t as the mean, i.e. $\mathbb{E}[A^d (t)] = \omega t$ and $t^{-1}\sqrt{\mathrm{Var}(A^d (t))} \rightarrow \sqrt{\beta \gamma}$ as $t \rightarrow \infty$. Thus, $\mathbb{P}(A^d (t) > 0)$ does not approach 0 or 1 as t increases.

5. Heavy-traffic limit for the $\sum_{i=1}^{n} P_i/GI/1$ queue

We now consider the single-server queue with arrival process $A^{d,n} (t) = A^{n} (\mu_n t)$ defined before (15). We assume that the service times are independent of the arrival process, mutually i.i.d. with a general distribution having mean $1/\beta$ and squared coefficient of variation (variance divided by the square of the mean) $c_s^2$, where the service times are independent of the arrival times. We work with the associated renewal counting process C(t). Since we center with $\beta n t$ in (15), we assume that the rate of this renewal process is also $\beta$. Let the scaled service process be

(17) \begin{equation}S_n (t) \equiv n^{-1/2}(C (nt) - \beta nt), \qquad t \ge 0,\end{equation}

As in [Reference Whitt36, Section 9.3], for the service process we only require a standard FCLT. Thus, that part of the following theorem can easily be generalized.

Let $X_n \equiv A^{d}_n - S_n$, $n \ge 1$, and let

(18) \begin{equation}Q_n (t) \equiv n^{-1/2}Q^n (t), \qquad t \ge 0,\end{equation}

where $\{Q^n (t)\,{:}\, t \ge 0\}$ is the queue length (number in system) process in the system with initial queue length $Q^n (0)$, arrival process $\{A^{d}_n (t)\,{:}\, t \ge 0\}$ defined before (15) and scaled service renewal counting process $\{S_{n} (t)\,{:}\, t \ge 0\}$ defined in (17).

With those definitions, we apply the one-dimensional reflection map in [Reference Whitt36, Section 13.5]. Let $\phi\,{:}\, D \times R \rightarrow D$ be the reflection map mapping a net input function x with $x(0) = 0$ and an initial queue length q(0) with $q(0) \ge 0$ into q(t) for $t > 0$ by

(19) \begin{equation}\phi (x) (t, q(0)) = q(0) + x(t) - \inf_{0 \le s \le t}{\{\min{\{q(0) + x(s), 0\}}\}}, \qquad t \ge 0.\end{equation}

The reflection map in (19) is a continuous function on its domain.

Proof. We apply the standard methodology for establishing a heavy-traffic FCLT for a single-server queue. We apply Donsker’s FCLT for the service times in [Reference Whitt36, Section 4.3] and the inverse equivalence in [Reference Whitt36, Theorem 7.3.2], in particular [Reference Whitt36, Corollary 7.3.2 on p. 236], to obtain $S_n \Rightarrow \beta^{3/2} c_s B$ in D as $n \rightarrow \infty$, where B is a standard Brownian motion. Then we can apply Corollary 3 to obtain the limit $A^{d}_n \Rightarrow A + \omega e$ in D. Joint convergence for $(A^{d}_n, S_n)$ then follows from independence and [Reference Whitt36, Theorem 11.4.4]. We can then apply [Reference Whitt36, Theorems 9.3.3, 9.3.4, 9.8.2].

By the continuous mapping theorem with the projection map, we have the following corollary providing a limit for the marginal distributions. Theorem 5 of [Reference Fendick16] provides the explicit form of the marginal distribution of the limit process, so that it can provide useful numerical results. Let the PDF of the joint limiting distribution be denoted by $f(x_{s}, q_{s},q_{s + t}) \equiv f_{X(s), Q(s), Q(s + t)} (x_{s}, q_{s},q_{s + t})$, and similarly for the associated marginal PDFs and conditional PDFs. We express the limiting distribution in terms of the exponential function and the standard normal CDF $\Phi (x) \equiv \mathbb{P}(N(0,1) \le x)$ and PDF $\phi (x)$. Let the associated CDF of $N(m, \sigma^2) \stackrel{\rm d}{=} m + \sigma N(0,1)$ be denoted by $\Phi (x;\ m, \sigma^2)$, and similarly for the others. To connect with [Reference Fendick16], let

(20) \begin{equation}\omega_s \equiv \frac{\alpha^* \omega - \beta^* x_s}{\alpha^* - \beta^* s} , \qquad \beta_{s}^* \equiv \frac{\alpha^* \beta^*}{\alpha^* - \beta^* s},\end{equation}

where

(21) \begin{equation}\alpha^* \equiv \beta , \qquad \beta^* \equiv - \beta \gamma ,\end{equation}

as in Remark 4. Let $\delta (\!\cdot\!)$ be the Dirac delta function and let $\textbf{1}_A$ be the indicator function, equal to 1 on A and 0 elsewhere.

Corollary 6. (Marginal limiting distributions.) Under the conditions of Theorem 5, $(X_n (s), Q_n (s), Q_n (s + t)) \Rightarrow (X (s), Q(s), Q(s+t))$ in ${\mathbb R}^3$ as $n \rightarrow \infty$, where X(s) is a (mean-$\omega s$, variance $s(\alpha^* - \beta_{s}^* s))$ Gaussian random variable for $\beta_{s}^*$ in (20), while the joint limiting distribution has joint PDF

(22) \begin{equation}f(x_{s}, q_{s},q_{s + t}) = f(x_{s}) f(q_{s} \mid x_{s}) f(q_{s + t} \mid x_{s}, q_{s}),\end{equation}

where, assuming $\mathbb{P}(Q(0) = q_0) = 1$,

\begin{align*} f(x_{s}) & \equiv \phi(x_s;\ \omega s, s(\alpha^* - \beta^* s)) = \frac{1}{\sqrt{s(\alpha^* - \beta^{*} s)}\,} \phi\left( \frac{x_s - \omega s}{\sqrt{s(\alpha^* - \beta^* s)}\,}\right)\!, \\[3pt] f(q_{s} \mid x_{s}) & \equiv \left(1 - \mathrm{e}^{\{-2q_{s}(q_{s} - x_s)/(\alpha^*s) \}}\right)(\delta(q_s - q_{0} - x_s)\textbf{1}_{\{q_s \ge 0\}} + \textbf{1}_{\{q_s - q_0 - x_s \ge 0\}}\delta(q_s)) \\[3pt] & \quad + \left(\frac{(4 q_s - 2 x_s)\mathrm{e}^{\{- 2q_s( q_s - x_s)/(\alpha^* s)\}}}{\alpha^* s}\right) \textbf{1}_{\{q_{s} - q_0 - x_s \ge 0\}}\textbf{1}_{\{q_s \ge 0\}}, \\[3pt] f(q_{s + t} \mid x_{s}, q_{s}) & \equiv \frac{1}{\sqrt{t(\alpha^* - \beta_{s}^{*} t)}\,}\phi\left( \frac{q_{s+t} - q_{s} - \omega_{s} t}{\sqrt{t(\alpha^* - \beta_{s}^* t)}\,}\right) \\[3pt] & \quad + \mathrm{e}^{\{-2 q_{s+t}(\beta_{s}^* q_{s + t} - \alpha^* \omega_s)/\alpha^{* 2}\}} (A_1 + A_2) \end{align*}

for

\begin{align*} A_1 & \equiv \left(\frac{4 \beta_{s}^* q_{s + t} - 2 \alpha^* \omega_s}{\alpha^{*2}}\right) \Phi\left( \frac{\left(2 \beta_{s}^* q_{s+t} - \alpha^* \omega_{s}\right)t - \alpha^*(q_{s+t} + q_{s})} {\alpha^* \sqrt{t(\alpha^* - \beta_{s}^* t)}}\right) , \\[3pt] A_2 & \equiv \left(\frac{\alpha^* - 2 \beta_{s}^* t}{\alpha^*\sqrt{t(\alpha^* - \beta_{s}^* t)}\,}\right)\phi\left( \frac{\left(2 \beta_{s}^* q_{s+t} - \alpha^* \omega_{s}\right)t - \alpha^*(q_{s+t} + q_{s})}{\alpha^* \sqrt{t(\alpha^* - \beta_{s}^* t)}}\right) , \end{align*}

$q_s, q_{s+t} \ge 0$. As a consequence, the associated conditional CDF, for $q_{s+t} \ge 0$, is

(23) \begin{align} \mathbb{P}(Q(s + t) \le q_{s+t} \mid X(s) &= x_s, Q(s) = q_{s}) = \Phi\left( \frac{q_{s+t} - q_{s} - \omega_{s} t}{\sqrt{t(\alpha^* - \beta_{s}^* t)}}\right)\nonumber\\[3pt] &\quad - \mathrm{e}^{\left(\frac{ -2 q_{s+t}(\beta_{s}^* q_{s + t} - \alpha^* \omega_s)}{\alpha^{* 2}}\right)}\Phi\left(\frac{(2 \beta_{s}^* q_{s+t} - \alpha^* \omega_{s})t- \alpha^*(q_{s+t} + q_{s})}{\alpha^*\sqrt{t(\alpha^* - \beta_{s}^* t)}}\right) \end{align}

and, for $q_t \ge 0$,

\begin{multline*} \mathbb{P}(Q(t) \le q_{t}) = \Phi\left( \frac{q_{t} - q_{0}- \omega t}{\sqrt{t(\alpha^* - \beta^* t)}\,}\right) \\[3pt] - \mathrm{e}^{\left(\frac{ -2 q_{t}(\beta^* q_{t} - \alpha^* \omega)}{\alpha^{* 2}}\right)}\Phi\left(\frac{(2 \beta^* q_{t} - \alpha^* \omega)t- \alpha^*(q_{t} + q_{0})}{\alpha^*\sqrt{t(\alpha^* - \beta^* t)}}\right)\!.\end{multline*}

Proof. We give a brief overview of the proof in [Reference Fendick16] (appearing on the last two pages), which draws heavily on [Reference Hajek22]. We focus on $X(t) \equiv Y(t) + \omega t$, where Y is an $(\alpha^*,\beta^*)$ $\Psi$-GMP for $(\alpha^*,\beta^*) = (\beta + \beta^3 c_s^2, - \beta \gamma)$ and $\omega t$ is the drift, as in Theorem 5. We exploit known results for the Brownian bridge by looking at increments from the past conditioned on later process values. For that purpose, let $X^{(s)} (t) \equiv X(s + t) - X(s)$ for some (s, t) with $0 \le s \le t$. We observe that conditioning $X^{(s)}$ on both $X(s) = x_s$ and $Q(s) = q_s$ results in a new $(\alpha^*,\beta_{s}^*)$ $\Psi$-GMP for $\beta_{s}^*$ in (20) that depends on $x_s$ but not $q_s$. Further conditioning it on $X^{(s)} (t) = x^{(s)}_t$ for some $t\ge s$ results in yet another $(\alpha^*, t^{-1} \alpha^*)$ $\Psi$-GMP with drift $t^{-1} x^{(s)}_t$ on [0, t] (which no longer depends on $\beta_{s}^*$). That process depends on $x^{(s)}_t$ but on neither $q_s$ nor $x_s$. Therefore, the process obtained by conditioning $\{Q(u)\,{:}\, s \le u \le s + t\}$ on $X(s) = x_s$, $Q(s) = q_s$ and $X^{(s)} (t) = x^{(s)}_t$ begins at state $q_s$ at time s and evolves according to a net input process that is an $(\alpha^*, t^{-1} \alpha^*)$ $\Psi$-GMP on the interval [0, t]. That process is a scaled Brownian bridge, for which the distribution of the queue length was previously obtained by [Reference Hajek22]. The result in (23) is obtained from that conditional queue length distribution using the law of total probability. We remark that $f(q_s \mid x_s)$ in (22) is the density of the CDF

\begin{displaymath} \mathbb{P}(Q (s) \le q_s \mid X (s) = x_s) = \left(1 - \mathrm{e}^{\{- 2q_s( q_s - x_s)/(\alpha^* s)\}}\right) \textbf{1}_{\{q_{s} - q_0 - x_s \ge 0\}}\textbf{1}_{\{q_s \ge 0\}} \end{displaymath}

as derived by [Reference Hajek22]. It suffices to apply the product rule for differentiation in this equation.

Some further insight can be gained from further Gaussian LLN limits for the limit process. For that purpose, let $X_s (t) \equiv (X(s + t) \mid X(s)= x_s)$ and $Q_s (t) \equiv (Q(s + t) \mid X(s) = x_s,Q(s) = q_s)$. Let $N(m, \sigma^2)$ denote a normal random variable with mean m and variance $\sigma^2$, and let $(x)^{+} \equiv \max{\{x, 0\}}$, so that $\mathbb{P}(N(m, \sigma^2)^{+} = 0) = \mathbb{P}(N(m, \sigma^2) \le 0)$.

Corollary 7. (LLN limit for the limit process.) Given the $\Psi$-GMP limit in Theorem 4 and the associated heavy-traffic limit in Theorem 5 and Corollary 6, we have the following limit:

(24) \begin{multline} t^{-1}(A^d (t), S(t), X(t), Q (t), X_s (t), Q_{s}(t)) \\ \Rightarrow (N_1(\omega, \beta \gamma), 0, N_1(\omega, \beta \gamma), N_1(\omega, \beta \gamma)^{+}, N_2(\omega_{s}, -\beta^{*}_{s}), N_2(\omega_{s}, -\beta^{*}_{s})^{+}) \text{ in } {\mathbb R}^6 \end{multline}

as $t \rightarrow \infty$ for some constant $s > 0$ and $(\omega_{s}, \beta^{*}_{s})$ in (21), where $(N_1, N_2)$ is a random vector in ${\mathbb R}^2$ with Gaussian one-dimensional marginal distributions.

Proof. For the first three processes, we exploit the distribution as a function of t as in the proof of Corollary 5. For the next two processes, note that $\mathbb{P}(X(t) \le xt) \rightarrow \Phi((x - \omega)/\sqrt{\beta \gamma})$ as $t \rightarrow \infty$ as just described, and

\begin{eqnarray*}\mathbb{P}(Q (t) \le t q \mid X(t) = tx) & = & (1 - \mathrm{e}^{-(2t^2q(q-x))/\alpha^* s})\textbf{1}_{\{tq - tx \ge q_0\}}\textbf{1}_{\{q \ge 0\}} \nonumber \\ & = & (1 - \mathrm{e}^{-(2t^2q(q-x))/\alpha^* s})\textbf{1}_{\{t \ge q_0/(q-x)\}}\textbf{1}_{\{q - x > 0\}}\textbf{1}_{\{q \ge 0\}} \nonumber \\ & \rightarrow & \textbf{1}_{\{q - x > 0\}}\textbf{1}_{\{q \ge 0\}} \qquad \text{ as } t \rightarrow \infty.\end{eqnarray*}

The joint PDF for the two limits is therefore

\begin{displaymath}\bar{f}(x,q) \equiv \bar{f}(x)\bar{f}(q \mid x) = \frac{1}{\sqrt{\beta \gamma}\,}\phi((x - \omega)/\sqrt{\beta \gamma})(\delta(q-x)\textbf{1}_{\{q \ge 0\}} + \delta(q)\textbf{1}_{\{q > x\}}).\end{displaymath}

The joint PDF for the limits of the last two processes is similarly obtained.

Figure 2 illustrates the heavy-traffic limits established in this section by plotting results of simulations of the scaled queue-length process in the $\sum_{i = 1}^{n} P_i/D/1$ queue.

Figure 2: Display of the empirical CDF of $Q_n (t)/t$ for $0 \le t \le T$ in the $\sum P_i/D/1$ queue with $Q_n (t)$ in (18), $\beta=\gamma=\mu = 1$ for $T = 25$ (left) and 500 (right), compared to $Q(t)/t$ in (24).

We consider the $(\beta,\gamma)$ $\Psi$-GPP with parameters $\beta = \gamma = 1$ as a function of the number n of processes superposed. We introduce a negative drift in the limit process by using (14) with $\mu = 1$. We let the service times initially have mean (and fixed value) 1, but that becomes $1/n$ with the scaling in (17). We look at n ranging from 4 to 200. We consider two values of T: 25 and 500. The similarity of the plots supports Corollary 7. The convergence as n increases supports Corollary 6. Numerical study indicates errors roughly of order $1/\sqrt{n}$.

6. Steady-state results for queues with $\Psi$-GPP arrivals

This section is something of a departure from Sections 3–5, because we no longer consider FCLTs or heavy-traffic limits. In this section we present further results for queues with $\Psi$-GPP arrivals. In Section 6.1 we obtain some elementary stability results and show that we can obtain both positive and negative results from a sample-path version of Little’s law. In Section 6.2 we show that stability can be achieved when there are infinitely many servers. In Section 6.3 we show that it is also possible to devise adaptive service policies that stabilize performance in a single-server queue with $\Psi$-GPP arrivals. We use new (but similar) notation in this section.

6.1. Stability and Little’s law

We first discuss the implications of the LLN in Theorem 3 for the $P/GI/1$ queue with a Pólya arrival process. For this purpose, let the service times have mean 1. Given Theorem 3, the following is a standard heavy-traffic law of large numbers for overloaded queues, as in [Reference Whitt36, Theorem 5.3.2].

Corollary 8. (Explosion.) Let Q(t) be the queue length process, starting with $Q(0) = 0$, in the $P/GI/1$ queue with a $(\beta, \gamma)$ $\Psi$-GPP arrival process, where $0 < \beta < \infty$ and $0 < \gamma < \infty$, and i.i.d. service times with mean 1. Then $t^{-1} Q(t) \rightarrow \max{\{L(\beta, \gamma) - 1, 0\}}$ as $t \rightarrow \infty$ with probability 1, so that $\mathbb{P}(Q(t) \rightarrow \infty \text{ as } t \rightarrow \infty) = \mathbb{P}(L(\gamma,\beta) >1)$, where $0 < \mathbb{P}(L(\gamma, \beta)\break > 1) < 1$.

We illustrate the less critical role of the usual traffic intensity by showing the results of a simulation experiment. Figure 3 plots the 25 individual sample paths of $Q(t)/t$ over [0, 250], starting empty, for a $P/D/1$ queue with parameter pairs $(\beta, \gamma) = (0.5, 1)$ (left) and $(\beta, \gamma) = (1.5, 1)$ (right). These plots only clearly display the paths that tend to be unstable by time $t = 250$. There are 3 of these in the 25 plots for $\beta = 0.5$ on the left, but about 12 in the 25 plots for $\beta = 1.5$ on the right. However, Corollary 2 implies that we cannot draw a firm conclusion about ultimate stability by what we see over the interval [0, 250].

Figure 3: Display of 25 individual sample paths of $Q(t)/t$ for $0 \le t \le 250$, starting empty, for a $P/D/1$ queue with parameter pairs $(\beta, \gamma) = (0.5, 1)$ (left) and $(\beta, \gamma) = (1.5, 1)$ (right).

Nevertheless, there is a version of Little’s law in this setting that can accommodate quite general initial conditions, but for simplicity assume that the system starts empty. Let $W_k$ be the waiting time of customer k, and let

(25) \begin{equation}L^* \equiv \lim_{t \rightarrow \infty}{t^{-1} \int_{0}^{t} Q(s) \, \mathrm{d} s}, \qquad W^* \equiv \lim_{n \rightarrow \infty}{n^{-1} \sum_{i = 1}^{n} W_k} , \qquad \lambda^* \equiv L(\beta,\gamma).\end{equation}

Corollary 9. (Little’s law.) Consider the $P/GI/1$ queue with a $(\beta, \gamma)$ $\Psi$-GPP arrival process starting with $Q(0) = 0$ as above, where $0 < \beta < \infty$ and $0 < \gamma < \infty$.

$(\textit{a})$ If $W^*$ in (25) is well defined and finite and if $L(\beta,\gamma) < 1$, which occurs with positive probability, then $L^*$ in (25) is well defined and $L^* = \lambda^* W^* < \infty$.

$(\textit{b})$ If $L(\beta, \gamma) > 1$, as occurs with positive probability, then $L^*$ in (25) is well defined but $L^* = \infty$. Moreover, $W^*$ is then also well defined and $W^* = L^* = \infty$.

Proof. For part (a), we apply the sample-path version of Little’s law from [Reference Stidham35]. For part (b), we apply Corollary 8 together with [Reference Wolff and Yao40, (6)]. Sufficient conditions for more results are also given in [Reference Wolff and Yao40].

6.2. The $\Psi$-GPP/GI/$\infty$ queue

We next consider the $\Psi$-GPP/GI/$\infty$ infinite-server queue. We show that bounded service times ensure reaching steady state in finite time. Let $\stackrel{\rm d}{=}$ denote equality in distribution.

Theorem 6. (Infinite-server queue with bounded service times.) For an infinite-server queue with a $(\beta, \gamma)$ $\Psi$-GPP arrival process, if the service times are (i) independent of the arrival process, (ii) mutually i.i.d., each distributed as a random variable V, and (iii) $\mathbb{P}(V \le \zeta) = 1$ for some $\zeta < \infty$, then the number of busy servers Q(t) reaches steady state by time $\zeta$, i.e.

(26) \begin{equation}Q(u) \stackrel{\rm d}{=} Q(t) , \quad \mathbb{E}[Q(u)] = \beta \mathbb{E}[V] \qquad \text{for all } u \ge t \ge \zeta.\end{equation}

Proof. We apply a useful device from [Reference Glynn and Whitt19, Section 2]. We note that if the service times are deterministic with $\mathbb{P}(V = d) = 1$, then

(27) \begin{equation}Q(t) = N(t) - N(t-d) \qquad \text{for all } t \ge d.\end{equation}

Because the $\Psi$-GPP arrival process has stationary increments, the distribution of Q(t) in (27) is independent of t for $t \ge d$, and so has reached steady state by time d.

We next observe that Q(t) also reaches steady state in finite time when the service time distribution has finite support within the interval $[0, \zeta]$. Hence, suppose that $P(V = d_i) = q_i$, $1 \le i \le n$. Now we can classify the arrivals by ‘type’ according to their service time, where this type assignment is done independently of the arrival process. Thus, the distribution of Q(t) is the sum of the number of each type in the system at time t, so that we can write

\begin{displaymath} Q (t) = \sum_{i = 1}^{n} Q_i (t) = \sum_{i = 1}^{n} [N(t) - N(t - d_i)] T_i (t, d_i) ,\end{displaymath}

where $T_i (t, d_i)$ is the proportion of the $N(t) - N(t - d_i)$ arrivals in $[t - d_i, t]$ that are of type i, i.e. have service times $d_i$, which has a multinomial distribution.

We then observe that the distribution of Q(t) is independent of t, because the random vector $(N(t) - N(t - d_i)\,{:}\, 1\le i \le n)$ is independent of t in ${\mathbb R}^n$ for all $t \ge \max{\{d_i\,{:}\, 1 \le i \le n\}}$.

The joint distribution of $(T_i (t, d_i)\,{:}\, 1 \le i \le n)$ is somewhat complicated, but we can directly deduce that it has property (26) and we can write down the mean:

\begin{eqnarray*} \mathbb{E}[Q (t)] & = & \sum_{i = 1}^{n} \mathbb{E}[Q_i (t)] = \sum_{i = 1}^{n} \mathbb{E}[N(t) - N(t - d_i)] \mathbb{E}[T_i (t, d_i)] \nonumber \\& = &\sum_{i = 1}^{n} \beta d_i \mathbb{E}[T_i (t, d_i)] = \beta \mathbb{E}[V].\end{eqnarray*}

Since the probability distributions with finite support in $[0, \zeta]$ are dense in the space of all probability distributions with support in $[0, \zeta]$, we obtain the general results by taking a limit. To obtain the limit, note that we can express Q(t) in terms of the arrival times $T_i$ associated with the arrival process N(t) and service times $V_i$ by $Q(t) = \sum_{i = 1}^{N(t)} \textbf{1}_{\{T_i + V_i > t\}}$, which representation implies that Q(t) is almost surely a continuous function of the service times with respect to the limit process, so that we can apply the generalized continuous mapping theorem from [Reference Whitt36, Theorem 3.4.4]. In particular, we see that Q(t) is almost surely a continuous function of the service times except when $T_i + V_i = t$ for some i. For any CDF of V, that almost surely does not occur because

\begin{displaymath}\sum_{i = 1}^{\infty} \mathbb{P}(t \le T_i \le t + \varepsilon) =\mathbb{P}(N(t+ \varepsilon) - N(t) \ge 1) \le \mathbb{E}[N(t+ \varepsilon) - N(t)] = \beta \varepsilon\end{displaymath}

by (3). We see that the probability is 0 by letting $\varepsilon \downarrow 0$. Hence we can apply the generalized continuous mapping theorem to deduce that the distribution of Q(t) is almost surely a continuous function of the distribution of V in this setting.

6.3. Adaptive service processes to enforce stability

Suppose that we let the service rate at time t depend on the history of the arrival process up to that time. In particular, let the arrival-history-dependent service rate be

(28) \begin{eqnarray}\mu (t) \equiv \mu(t \mid {\mathcal{H}}_t) = \frac{\lambda (t)}{\rho} = \frac{\gamma N(t) + \beta}{\rho(\gamma (t) + 1)} \qquad \text{for all } t \ge 0.\end{eqnarray}

The proposed control is a variant of the rate-matching service-rate control in [Reference Ma and Whitt30, Reference Whitt37]. The following is an analog of [Reference Whitt37, Theorem 3.1].

Theorem 7. If the arrival-history-dependent service process in (28) with $\rho < 1$ is used in the single-server queue with $Q(0) = 0$ and the GPP arrival process having parameter triple $(\kappa (t), \gamma, \beta)$ with $\kappa$ in (2), then the pair of rates satisfies the joint limit

(29) \begin{equation}(\lambda (t), \mu (t)) \rightarrow (L(\gamma, \beta), L(\gamma, \beta)/\rho) \qquad \text{as } t \rightarrow \infty\ \textit{with\ probability\ 1} , \end{equation}

where $L(\gamma, \beta)$ is the gamma random variable in Theorem 3. Moreover, the queue length process $\{Q(t)\,{:}\, t \ge 0\}$ is distributed the same as in an $M/M/1$ queue in a random environment, so that

(30) \begin{equation}Q(t) \Rightarrow Q \text{ as } t \rightarrow \infty, \qquad {where \,} \mathbb{P}(Q = k) = (1 - \rho) \rho^k, \quad k \ge 0,\end{equation}

as in the $M/M/1$ queue with traffic intensity (arrival rate divided by the maximum potential service rate) $\rho$.

Proof. First, the joint limit in (29) follows easily from Corollary 1 and (28). For the connection to $M/M/1$, we start by representing the counting processes N(t) and S(t) in terms of Poisson processes, using the martingale representation for point processes as in [Reference Baccelli and Bremaud3, Section 1.8]. Let $\Pi_1 (t)$ and $\Pi_{2} (t)$ be two independent rate-1 Poisson counting processes. Then the counting process N can be represented as $\Pi_1 (\Lambda (t))$, $t \ge 0$, where $\Lambda (t)$ is the cumulative rate function associated with $\lambda (t)$ and the compensator of the point process, i.e. $\Lambda (t) \equiv \int_{0}^{t} \lambda (s) \, \mathrm{d} s$, $t \ge 0$. Because of (28), we thus have $\{(N(t), S(t))\,{:}\, t \ge 0\} \stackrel{\rm d}{=} \{(\Pi_{1} (\Lambda (t)), \Pi_{2} (\rho^{-1} \Lambda (t)))\,{:}\, t \ge 0\}$; see, e.g., [Reference Ethier and Kurtz13, Theorem 6.4.1]. Then the queue length process Q(t) can be expressed as the reflection map applied to the net-input process $Y(t) \equiv N(t) - S(t)$. Let $\tilde{Q} (t)$ and $\tilde{Y} (t)$ be the corresponding processes in an $M/M/1$ queue with arrival rate 1 and service rate $1/\rho$, i.e. $\tilde{Y} (t) \equiv \Pi_{1} (t) - \Pi_{2} (\rho^{-1} t)$. With this representation we have $\{Y(t), Q(t)\,{:}\, t \ge 0\} \stackrel{\rm d}{=} \{(\tilde{Y} (\Lambda (t)), \tilde{Q} (\Lambda (t)))\,{:}\, t \ge 0\}$, where $\{\Lambda (t)\,{:}\, t \ge 0\}$ is independent of $\{(\tilde{Y} (t), \tilde{Q} (t))\,{:}\, t \ge 0\}$. Because $\Lambda (t) \rightarrow \infty$ as $t \rightarrow \infty$ with probability 1, we have the claimed representation and convergence in (30).

7. Completing the proof of Theorem 4

We complete the proof of Theorem 4 by verifying the two inequality conditions in Hahn’s theorem [Reference Hahn21] as in [Reference Whitt36, Theorem 7.2.1]. To do so, we prove convergence in D[0, 1], as in [Reference Hahn21], but we note that essentially the same arguments show convergence in D[0, T] for any $T>0$ and therefore implies convergence in $D[0,\infty)$; see [Reference Whitt36, Section 12.9].

By Theorem 1, $\mathbb{E}[A_1(s)A_1(t)]= \beta s(1+\gamma t)$ when $A_1$ is the centered process defined in (9). First, for $0\le t_1\le t_2\le 1$,

\begin{align*} \mathbb{E}\left[(A_1(t_2)-A_1(t_1))^2\right] &=\beta \left(t_2-t_1\right)\left(1+\gamma \left(t_2-t_1\right)\right) \nonumber \\& = \beta \left(t_2+\gamma t_2^2\right)-\beta \left(t_1+\gamma t_1^2\right)-2\beta \gamma t_1\left(t_2-t_1\right)\nonumber \\& \le \beta \left(t_2+\gamma t_2^2\right)-\beta \left(t_1+\gamma t_1^2\right) , \end{align*}

so that condition (2.3) of [Reference Whitt36, Theorem 7.2.1] is met. We will also show that

(31) \begin{equation}\mathbb{E}[\left(A_1(t)-A_1\left(t_1\right)\right)^2\left(A_1\left(t_2\right)-A_1(t)\right)^2] \le c\left(t_2-t_1\right)^2\hfill\end{equation}

for $0\le t_1\le t\le t_2\le 1$ and a constant c, so that condition (2.4) of [Reference Whitt36, Theorem 7.2.1] will be met as well.

To do the explicit calculations, we again apply [Reference Cha7, Theorem 3 and Remark 3], as in the proof of Theorem 1. That applies for any GPP $N^1$, the conditional distribution for the sequence of times when $N^1$ increases on $\left(0,1\right)$ given that $N^1\left(1\right)=k$ is the same as that of the order statistics of k i.i.d. random variables, each with probability density function

\begin{equation*}f(x)\equiv \frac{\gamma \kappa (x)\exp\!\left(\gamma K(x)\right)}{\exp\!\left(\gamma K\left(1\right)\right)-1}, \qquad 0\le x\le 1.\end{equation*}

For a $\Psi$-GPP, $f(x)=1$. Hence, the conditional distribution for the ordered sequence of times when the $\Psi$-GPP $N^1$ increases on $\left(0,1\right)$ given that $N^1\left(1\right)=k$ is the same as that of the order statistics of i.i.d. uniform random variables $U_j$ on [0, 1].

Following the proof of [Reference Billingsley5, Theorem 14.3], let $p_1=t-t_1$, $p_2=t-t_2$; let $V_j$ be $\left(1-p_1\right)$ or $- p_1$ as $U_j$ lies in $\left[t_1,t\right)$ or not; let $W_j$ be $\left(1-p_2\right)$ or $- p_2$ as $U_j$ lies in $\left[t,t_2\right]$ or not. Also, let $V=\overset k{\underset{j=1}{\sum }}V_j$, $W=\overset k{\underset{j=1}{\sum }}W_j$, $\tilde{V}= \left(N^{1}(1) - \beta\right)p_1$ and $\tilde{W} =\left(N^{1}(1) - \beta\right)p_2$. Then

(32) \begin{align}&\mathbb{E}\left[\left(A_1(t)-A_1\left(t_1\right)\right)^2\left(A_1\left(t_2\right)-A_1(t)\right)^2 \mid N^1\left(1\right)=k\right] \nonumber \\[2pt] & \quad= \mathbb{E}\left[\left(\overset k{\underset{j=1}{\sum }}\left(V_j+\frac{\left(k-\beta\right)p_1} k\right)\right)^2\left(\overset k{\underset{j=1}{\sum}}\left(W_j+\frac{ \left(k-\beta\right)p_2} k\right)\right)^2 \Bigg| N^1\left(1\right)=k\right] \nonumber \\[2pt] &\quad = \mathbb{E}\left[(V + \tilde{V})^2(W + \tilde{W})^2 \mid N^1(1) = k\right] \nonumber \\[2pt]&\quad = \mathbb{E}\left[\left((VW + V\tilde{W}) +W\tilde{V} + \tilde{V}\tilde{W}\right)^2 \mid N^1(1) = k\right] \nonumber \\[2pt]&\quad \le 4 \mathbb{E}\left[(VW)^2 + (V\tilde{W})^2 + (W\tilde{V})^2 + (\tilde{V}\tilde{W})^2 \mid N^1(1) = k\right],\end{align}

where the last inequality follows from the Cauchy–Schwartz inequality.

It follows from [Reference Cha7, Theorem 1(i)] that $N\equiv N^1\left(1\right)$ has a negative binomial distribution with moment-generating function

\begin{equation*}M\left(s\right)=\left(\frac{\theta }{1-\left(1-\theta \right)\mathrm{e}^s}\right)^r ,\end{equation*}

where $r=\beta /\gamma $ and $\theta =\left(1+\gamma \right)^{-1}$. Then,

\begin{equation*}\mathbb{E}\left(N^p\right)=\left(\frac{\mathrm{d}^p}{\mathrm{d} s^p}M\left(s\right)\right)_{s=0} \qquad \text{for } p\ge 1 ,\end{equation*}

which implies that

(33) \begin{align} \begin{split} \mathbb{E}[N] & = \beta, \\[2pt] \mathbb{E}\left[N^2\right] & = \beta \gamma +\beta ^2+\beta , \\[2pt] \mathbb{E}\left[N^3\right] & = \beta (2\gamma ^2+3\beta \gamma +3\gamma +\beta ^2+3\beta +1) , \\[2pt] \mathbb{E}\left[N^4\right] & = \beta (6\gamma ^3+11\beta \gamma ^2+12\gamma ^2+6\beta ^2\gamma +18\beta \gamma +7\gamma +\beta^3+6\beta ^2+7\beta +1). \end{split}\end{align}

By [Reference Billingsley5, (14.10)], $\mathbb{E}\big[\left(VW\right)^2 \mid N^1\left(1\right)=k\big]\le 6\ k^2p_1p_2$, so that

(34) \begin{equation}\mathbb{E}\left[\left(VW\right)^2\right]\le 6 E\left[N^2\right]p_1p_2.\end{equation}

Similarly,

\begin{align*}\mathbb{E}\left[\left(V\tilde{W}\right)^2\right] & = \mathbb{E}\left[Np_1\left(1-p_1\right)\left(\left(N-\beta\right)p_2\right)^2\right], \\[2pt] \mathbb{E}\left[\left(\tilde{V} W\right)^2\right] & = \mathbb{E}\left[Np_2\left(1-p_2\right)\left(\left(N-\beta\right)p_1\right)^2\right],\end{align*}

and

(35) \begin{equation}\mathbb{E}\left[\left(\tilde{V}\tilde{W}\right)^2\right]= \mathbb{E}\left[\left(\left(N-\beta\right)p_1\right)^2\left(\left(N-\beta\right)p_2\right)^2\right].\end{equation}

Using Macsyma for algebraic simplification, we find first, by (32),

\begin{align*} &\mathbb{E}[\left(A_1(t)-A_1\left(t_1\right)\right)^2 \left(A_1\left(t_2\right)-A_1(t)\right)^2] \\[2pt] &\quad \le \Omega \equiv 4\left(\mathbb{E}\left[\left(VW\right)^2\right]+\mathbb{E}\left[\left(V\tilde{W}\right)^2\right]+\mathbb{E}\left[\left(\tilde{V} W\right)^2\right]+ \mathbb{E}\left[\left(\tilde{V}\tilde{W}\right)^2\right]\right)\!.\end{align*}

Then, by (33)–(35), we find that

\begin{align*} \Omega & = 4\beta\left(6 \gamma^3 p_1p_2 + 3\beta\gamma^2 p_1 p_2 + 8\gamma^2 p_1p_2 +4\beta \gamma p_1p_2+ \gamma p_1p_2 +\beta p_1p_2 - p_1p_2 + 2 \gamma^2 p_2 \right. \\[2pt] & \left. \quad + \, \beta \gamma p_2 + 3\gamma p_2+\beta p_2+ p_2 +2\gamma^2 p_1 + \beta \gamma p_1 + 3 \gamma p_1 + \beta p_1 + p_1 + 6\gamma + 6\beta + 6\right)p_1 p_2 \\[2pt] & \le 4\beta (6 \gamma^3 +3\beta\gamma^2 + 12 \gamma^2 + 6\beta\gamma +13 \gamma + 9\beta + 7)p_1p_2 \\[2pt] & = 4\beta (6 \gamma^3 +3\beta\gamma^2 + 12 \gamma^2 + 6\beta\gamma +13 \gamma + 9\beta + 7)\left(t-t_1\right)\left(t_2-t\right) \\[2pt] & \le 4\beta (6 \gamma^3 +3\beta\gamma^2 + 12 \gamma^2 + 6\beta\gamma +13 \gamma + 9\beta + 7)\left(t_2-t_1\right)^2.\end{align*}

Thus, we conclude that condition (2.4) of [Reference Whitt36, Theorem 7.2.1] holds for $c=4\beta (6 \gamma^3 +3\beta\gamma^2 + 12 \gamma^2 + 6\beta\gamma +13 \gamma + 9\beta + 7)$. Under those conditions, [Reference Whitt36, Theorem 7.2.1] shows that $A_n\Rightarrow A$, where A is a zero-mean Gaussian process with the same covariance kernel as $A_1$.