Queueing and risk models with dependencies

This paper analyzes various stochastic recursions that arise in queueing and insurance risk models with a ‘semi-linear’ dependence structure. For example, an interarrival time depends on the workload, or the capital, immediately after the previous arrival; or the service time of a customer depends on her waiting time. In each case, we derive and solve a fixed-point equation for the Laplace–Stieltjes transform of a key performance measure of the model, like waiting time or ruin time.


Introduction
One of the most fundamental relations in queueing theory is the Lindley recursion, which is a recursive relation between the waiting times W i and W i+1 of the ith and (i + 1)st customers in the single-server queue. With (B i ) i a sequence of service times, and (A i ) i a sequence of interarrival times, the sequence (W i ) i satisfies This paper is dedicated to Masakiyo Miyazawa in friendship and in admiration for his numerous deep contributions to queueing theory.
One alternatively says that (W i ) i is a random walk with increments (B i −A i ) i reflected at zero. The Lindley recursion has been studied in much generality when it comes to the specific distributional assumptions imposed on the random sequences (B i ) i and (A i ) i . However, in the vast majority of all papers it is assumed that both sequences consist of independent and identically distributed (i.i.d.) random variables, that both sequences are independent, and in addition, that all A i and B i are independent of all W j , for j ≤ i. The goal of the present paper is to explore a class of stochastic recursions (1) in which some of the above-mentioned independence assumptions are lifted and for which, nevertheless, a detailed exact analysis can be provided. Below we briefly discuss the models under consideration.
• In Sects. 2 and 3, we assume that A i = c(W i + B i ) + A i , where (A i ) i is an independent sequence of exponentially distributed random variables (in both sections some relaxation of the exponentiality assumption is possible) and where 0 < c < 1. We thus model a positive correlation between a customer's sojourn time W i + B i and the time until the next customer arrives. A possible application is a queueing system with a controller: If the workload becomes high after an arrival, the controller scales down the arrival rate of jobs. Alternatively, this scaling down could reflect the strategy of the customers; they might be reluctant to join the queue when the workload is high. The key performance measure under consideration in Sect. 2 is the steady-state waiting time; its LST (Laplace-Stieltjes transform) is presented in Theorem 2.1. We also derive an asymptotic expansion of the probability of zero waiting time, and of the waiting time moments, in the parameter c = ↓ 0 that indicates how close the model is to the classical M/G/1 queue. The key performance measure in Sect. 3 is the time τ (x) until the system becomes empty, when the initial capital is x. In this section, we also reinterpret the model as a so-called dual risk model. The dual risk model is used to study the possible ruin of companies, and accordingly, the time until the queue becomes empty translates into the time a company gets ruined when its initial capital is x. The Laplace transform (with respect to that initial capital x) of the ruin time LST Ee −sτ (x) , is derived in Theorem 3.1. • Section 4 studies the Cramér-Lundberg (CL) insurance risk model. In this model, generally distributed claims arrive according to a Poisson process, and in between claims, the capital of the insurance company increases at a fixed premium rate. The field of insurance risk has much in common with queueing theory; in particular, the CL model is known to be dual to the M/G/1 queue (cf. Chapter III.2 of [5]). We consider the following dependency structure in the CL model. The interarrival time between the ith and (i + 1)st claim depends on the capital y of the insurance company right after the ith claim, in very much the same way as the interarrival time depends on the workload in Sects. 2 and 3 (but taking into account that now jumps are downward and the linear slope is upward): When the capital, right after the ith claim, equals y, then the next interclaim time where (A i ) i is an independent sequence of exp(λ) distributed random variables. In Theorem 4.1, we determine the Laplace transform, with respect to initial capital x, of the ruin time LST Ee −sτ (x) .
• Finally, in Sect. 5, we turn to the single-server queue in which the service time of a customer depends, in a somewhat similar way as above, on her waiting time. We assume that B i = max(B i − cW i , 0), where (B i ) i is an independent sequence of exp(μ) distributed random variables. The negative correlation between waiting and service time includes the feature of zero service time if the waiting time is too large. We derive the steady-state waiting time and workload LST (Theorem 5.1).
All the above-sketched models exhibit a 'semi-linear' dependence structure, and in all these models, we arrive at a fixed-point equation for the LST ω(s) of the steady-state waiting time W , or for the Laplace transform of the ruin time LST, of the following form:  Hadidi [20] shows that the waiting times are hyperexponentially distributed, while Hadidi [21] studies the sensitivity of the waiting time distribution to the value of the correlation coefficient and Langaris [24] considers the busy period distribution. The case in which B i depends on the previous A i−1 has only been studied in a few cases; see the queueing studies [8,13,14]. For a broad discussion of dependence phenomena in queueing models of packet networks (dependence between successive interarrival times, or between successive service times, or between interarrival and service times), we refer to Fendick et al. [18]. A single-server queue in which the service rate of a customer with exponential service time distribution depends on her waiting time has been analyzed in [26]. There it is assumed that the service rate is μ(w) if the waiting time equals w, and assumes moreover that μ(w) is a step function. In insurance risk, more than in queueing, attention has been given to models with various forms of dependence. We refer to Chapter XIII of Asmussen and Albrecher [5] for an overview and to [3] (generalizing [23]) for a rather general class of Markovian risk models, in which claim interarrival times and claim sizes have a joint distribution that depends on the state of some underlying finite Markov chain. A similar queueing model with a Markov dependence structure was studied by Adan and Kulkarni [2]; see also [14]. A special case of the dependence structure in [3] is an insurance risk model in which an interarrival time depends on the previous claim size via a threshold mechanism; the queueing dual of that model was studied in [12]. We finally remark that, if one is satisfied with asymptotic results, more general forms of dependence may be allowed; see, e.g., [6,22].

The M/G/1 queue with interarrival times dependent on workload
In this section, we study the M/G/1 queue with the special feature that the time between the arrivals of customers n and n + 1 depends on the workload found by customer n. More specifically, if the workload just after the arrival of customer n equals x ≥ 0, then the next interarrival time equals cx plus an independent, exponentially distributed time. A brief model description (Sect. 2.1) is followed by an analysis of the LST of the steady-state workload (Sect. 2.2). Section 2.3 presents an asymptotic analysis for the case c ↓ 0; of course, in the limiting case of c = 0 the model reduces to the classical M/G/1 queue.

Model description
Consider the following variant of the classical M/G/1 queue. If the workload just after the ith arrival (i.e., the sojourn time) equals x ≥ 0, then the next interarrival time equals cx (for some c ∈ (0, 1); if c ≥ 1, then all waiting times are zero) increased by an independent, exponentially distributed term A i with mean 1/λ. This construction has the interesting feature that there is a positive correlation between a customer's sojourn time and the time until the next customer arrives. The successive service times constitute a sequence of i.i.d. random variables (B i ) i , which are also independent of the sequence (A i ) i ; their LST is denoted by β(·).

Analysis
Let W i denote the waiting time of the ith arriving customer. It is readily verified that the sequence (W i ) i obeys the recursion x + denoting max(0, x) (similarly, below x − denotes min(0, x)). In the following, we assume that the steady-state waiting time distribution exists, which we claim to hold for any c > 0 as long as EB 1 < ∞ (and actually even if E log(1 + B 1 ) < ∞; see the account of this issue for a very similar model in [10]). Now we let i → ∞ in (3) and study the limiting random variable W ; let ω(s) denote the LST corresponding to W . With A, B denoting generic random variables with the same distribution as A 1 , B 1 , we have, using that x + + x − = x and that A ∼ exp(λ), This fixed-point equation is very similar to the ones which were solved in [9,11]. The solution procedure is to iterate (4), successively expressing ω((1 − c) j s) into ω((1 − c) j+1 s), and to prove the convergence of the resulting infinite sum of products. It should be noticed that the motivation here is different from that in [9,11]: Here, we wish to model a dependence between sojourn time and next interarrival time. However, we can translate (3) into the equation In this way, by adapting Theorem 2.7 of [11], we obtain the following result. Define where Remark 2.2 It should be observed that one way to obtain the identity (6) is to deduce (5) and using (6).

Remark 2.4 A generalization of (3) is to let the ith interarrival time be given by
for some δ(·); put differently, we are considering a Lévy process. Equation (4) is now replaced by Again, there is a relation to an earlier study, namely [9], where one has to take service times with LST β(δ(s)) instead of β(s). (3) and considering the transient behavior of W n via its transform ∞ n=1 r n E[e −sW n |W 1 = w], the same recursion procedure will enable us to derive an expression for the latter transform; cf. Section 2.1 of [11] for a similar approach. We also refer to [11] for a generalization to the case in which the A i have a hypo-exponential distribution.

Asymptotic expansions
In this subsection, we are interested in the performance measures P(W = 0) and EW k , k = 1, 2, . . . in the regime that c is small, i.e., a perturbation of the classical M/G/1 queue. For such values of c, there is a rather weak dependence between sojourn time (workload just after an arrival) and the subsequent interarrival time, while c = 0 concerns the case without dependence. In view of the rather complicated form of the exact expressions for the above-mentioned performance measures, we are interested in obtaining good asymptotic expansions of them in the parameter c.

Remark 2.6
Our interest in these asymptotic expansions was strengthened by the fact that the French hydrologist E. Mouche has indicated (personal communication; cf. also [25]) that the model and analysis of [11] are relevant in certain hydrology models, in particular when a is close to one-which corresponds to our c being close to zero. It should be observed, though, that analyzing the model of [11] for a close to one is not precisely the same as analyzing the model of the present section for c close to zero, because the model translation also involves replacing B i by (1 − c)B i .
Our starting point for the asymptotic analysis for c = ↓ 0 is (4), which can be rewritten as Denote the steady-state waiting time when = 0, corresponding to the classical M/G/1 queue, byW . Assuming that the first K moments of W are well defined, we can give a Taylor series development of P(W = 0) and of those K moments of W up to m terms for some m ∈ N. In other words, for ↓ 0 and k ∈ {1, . . . , K }, Differentiating both sides of (8) with respect to s, substituting s = 0 and multiplying by −1, we obtain that and hence, with ρ := λ EB: Substituting (9), and (10) for k = 1, into (12) gives us R 0,1 : recognize in the latter expression the mean queue length in the M/G/1 queue. Furthermore, by equating the h powers in both sides of (12) for h = 2, 3, . . . , we obtain: Both sides of (4) are subsequently differentiated k = 2, . . . , K times with respect to s, after which s = 0 is substituted. After some calculations, we obtain, for k = 2, . . . , K , In particular, we recover the following familiar M/G/1 result [19] by inserting = 0: observe that from this expression all moments EW k can be computed recursively. Substituting (10) into (15) and using (16) gives, for k = 2, . . . , K : and hence, We return to our objective of identifying the coefficients R k,h .
• Equating factors on both sides allows us to express R k−1,1 into R k−2,1 , . . . , R 0,1 and known terms; regarding the latter, recall that all moments ofW can be obtained recursively from (16). We have, for k = 2, . . . , K : Since R 0,1 has been derived in (13), all R k,1 can be obtained. • By equating h factors, h = 2, 3, . . . , on both sides of (17), we can express R k−1,h into R k,h−1 and in terms R j,l with j + l ≤ k − 2 + h. We give the recursion for 2 terms: Considering the matrix (R g,i ) g≥0,i≥1 , we thus first obtain the element R 0,1 in the upper left corner; then the first column of elements R k,1 ; and then successively all the 'anti-diagonals' (where an anti-diagonal corresponds to the entries of which the sum of the indices is equal). In particular, R 1,1 gives R 0,2 , subsequently we use R 2,1 to obtain R 1,2 and R 0,3 ; then, we use R 3,1 to obtain R 2,2 , R 1,3 and R 0,4 ; etc.

The dual risk model
The classical Cramér-Lundberg insurance risk model is often viewed as being dual to the M/G/1 queue: The reserve level process grows linearly with the premium rate in between downward Poisson jumps, with these jumps (to be thought of as claim sizes) stemming from a general distribution. When there is equality between the arrival rates of the two Poisson processes, the two jump size distributions and the slopes, we have the remarkable property (cf. Section III.2 of [5]) that the ruin probability in the Cramér-Lundberg model when starting with initial capital x equals the probability that the steady-state waiting time in the M/G/1 queue exceeds x. However, in the insurance risk literature there is also some interest in a model, called the dual risk model, that has exactly the same sample path behavior as the M/G/1 queue, until level zero is first reached: jumps upward and a constant slope downward between jumps. Such a dual risk model represents the surplus of a company with a fixed expense rate and occasional gains; examples are pharmaceutical, R & D and petroleum companies. Various performance measures of the dual risk model have been studied since the seminal work by Avanzi et al. [7].
In the present section, we consider a dual risk model in which we make the exact same assumptions as in Sect. 2. It is evident that from a ruin probability point of view, nothing interesting happens: the ruin probability is 1 (see also Remark 3.2). However, it is interesting to determine the distribution of the time to ruin τ (x), when starting (immediately after an upward jump) at level x. Observe that this is equivalent with determining the distribution of the busy period of the model's M/G/1 counterpart, when starting at level x right after a jump (i.e., a customer arrival). Denote the LST of this ruin time τ (x) by K (s, x). Let T s denote an exponentially distributed random variable with mean 1/s. Observe that K (s, x) can be interpreted as the probability that the ruin time τ (x) occurs before T s .
Let us first restrict ourselves to the case of exp(μ) jumps upward, i.e., the corresponding distribution function is given by B(z) = 1 − e −μz . By looking ahead to the next jump, we can write  K (s, x) dx, we obtain by interchanging the integrations over x and y and substantial calculus that Writing Iteration of (22), which again has the form (2), finally results in the following theorem.
where an empty product is defined to be one and where the remaining unknown k(s, μ) is determined by substituting α = μ in (27).
In relation to the last statement of the theorem, observe that k(s, μ) appears with a prefactor in all G(s, ζ ( j) (α)) terms in (27). It is easily seen that the sum of products in (27) converges, because both G(s, ·) and H (s, ·) decrease geometrically fast (with a factor 1 − c) to zero. Finally, we would like to remark that μ k(s, μ) can be interpreted as the LST of a busy period starting from an empty system with an exp(μ) jump upward-or, equivalently, the LST of the ruin time in the dual risk model with exp(μ) distributed initial capital. K (0, x) equals the probability that ruin ever occurs, starting at level x. It is intuitively clear that this probability should be one for all c > 0. Indeed, substituting k(0, α) = 1/α and k(0, μ) = 1/μ in (22) gives an identity.

Remark 3.2 Notice that
It is readily verified that, for c = 0, (22) becomes Taking α = μ gives an identity, but observing that the term between square brackets in the left-hand side of (29) has exactly one zero in the right half α-plane, viz.
and that α 1 should also be a zero of the right-hand side of (29), quickly leads to (28).

Remark 3.4 See [16]
for a numerical recipe to numerically invert two-dimensional Laplace transforms; cf. also the practical guidelines presented in [4,17]. Using these techniques, one can reliably determine the density or distribution function of τ (x) from our expression for k(s, α).
So far in this section, we have assumed that the upward jumps are exponentially distributed. In the remainder of this section, we generalize this to the hyperexponential and Erlang cases.
• The hyperexponential case. In this case, the cumulative distribution function reads . It is quite straightforward to see that in this case The iterative approach that we applied for the exponential case in this section works in exactly the same way for the hyperexponential case. Formula (22) still holds, with now The only difference is that we now end up with K unknowns, namely k(s, μ i ), i = 1, . . . , K . These can be determined by substituting α = μ i , i = 1, . . . , K in the equivalent of (27).
• The Erlang case. In the case of an Erlang-K distribution, the cumulative distribution function reads In this Erlang-K case, too, it is straightforward to find k(s, α). Observing that we obtain The iterative approach that we applied for the exponential case can be applied in the exact same way for the Erlang-K case. Again Formula (22) still holds, with now The only difference is that the K unknowns we end up with are k(s, μ) and its first K − 1 derivatives. These are obtained by differentiating the equivalent of (27) 0, 1, . . . , K −1 times with respect to α, each time followed by substituting α = μ. This results in K linear equations for these K unknowns.

Remark 3.5
It is straightforward to combine these two cases further to the case of the distribution function B(·) corresponding to a weighted sum of Erlang-k distributions. It is well known that distributions from this class can approximate the distribution of any nonnegative random variable arbitrarily close.

The Cramér-Lundberg model with dependence between capital and interclaim times
In this section, we consider the following variant of the Cramér-Lundberg risk model. Claim sizes are i.i.d. with a non-specified cumulative distribution function B(·) and LST β(·), independent of everything else. In between claims, the capital grows at a constant rate, which is assumed (without losing any generality) to be one. When the capital right after the ith claim takes some value y ≥ 0, then the next interclaim time equals max(0, A i − cy), where c > 0 and (A i ) i are i.i.d. exp(λ) distributed random variables, independent of everything else.
It is directly seen that the underlying mechanism is such that a large capital gives rise to a relatively small interclaim time. This means that eventual ruin is certain in this model. (When the capital is large, there is a cascade of claims, so the capital is pulled down, whereas if the capital is small, the model effectively behaves as the conventional Cramér-Lundberg model.) Below we shall study the ruin time LST R(s, x) := E[e −sτ (x) ], where τ (x) is the ruin time when starting with capital x > 0 immediately after a claim.
Conditioning on the value of A i , and in particular distinguishing between this quantity being smaller than cx (leading to an immediate claim arrival) or larger than cx, we obtain, where Recognizing convolutions in both A I and A II , and remembering that the Laplace transform ofB(x) equals (1 − β(α))/α, we obtain and, after some interchange of integrations, with f (s, α) = s + λ(1 − c) + α, Notice that the zero α = s + λ(1 − c) of f (s, −α) is a removable singularity of A II , as it also makes the numerator of A II zero. It is readily verified that r (0, α) = 1/α, as it should. Combination of (38), (40) and (41) with M(s, α) : Iteration of (42) finally results in the following theorem.
It is easily seen that the above sum of products converges. Both N (s, α + iλc) and M(s, α + iλc) converge to zero as i → ∞; the convergence of the sum of products is at least geometrically fast for any α ≥ 0, as the absolute value of the product is bounded by a constant times [β(λc + α)] j .

The single-server queue with service time dependent on waiting time
This last section considers the following variant of the classical M/M/1 queue. Customers arrive according to a Poisson process with rate λ. If the waiting time W i of the ith arriving customer equals x ≥ 0, then her service time equals where c > 0 and where (B i ) i is a sequence of independent, exponentially distributed random variables with mean 1/μ, independent of anything else. Notice that there is a negative correlation between a customer's waiting time and her service requirement. In particular, when the waiting time is very large, the service requirement is likely to be zero. The probability of the latter event is P(B i < cW i ); if a customer has a positive service time, then it equals an exp(μ) quantity, due to the memoryless property of the exponential distribution.
an empty product denoting one. Then note that π 0 follows from φ(0) = 1: We conclude the following.

Theorem 5.1 The steady-state waiting time and workload LST is given by
(51) Apparently the distribution of W is an infinite weighted sum of exponential distributions with rates μ(1 + ic), i = 0, 1, . . . . Differentiating (51) and substituting s = 0 give the mean waiting time: When c = 0, φ(s) reduces to the familiar M/M/1 result where one has to require that λ < μ and where π 0 = 1 − λ/μ.

Remark 5.2
The transient behavior corresponding to the (W i ) i sequence can be found by following a similar procedure as in Section 2.1 of [11], cf. also Remark 2.5, yielding ∞ n=1 r n E[e −sW n | W 1 = w]. Notice that, after multiplication by (1 − r )/r , this also gives the LST of W N with N being geom(r ) distributed.

Remark 5.3
When the B i have a general distribution, the calculations in (47) become much more involved, and so far we have not been able to arrive at any tractable recursion. The hyperexponential case seems doable, though. If P(B > x) = K i=1 p i e −μ i x , then (47) should be replaced by and (48) generalizes to In [1], recursions of exactly this type are treated. The commutativity of ζ i (s) := s+μ i c and ζ j (s) := s + μ j c, i.e., ζ i (ζ j (s)) = ζ j (ζ i (s)), makes the recursion (54) relatively easy even though in each iteration step a term φ(s + K i=1 b i μ i c) gives rise to K new terms.