Functional equations with multiple recursive terms

In this paper, we study a functional equation for generating functions of the form f(z)=g(z)∑i=1Mpif(αi(z))+K(z)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f(z) = g(z) \sum _{i=1}^M p_i f(\alpha _i(z)) + K(z)$$\end{document}, viz. a recursion with multiple recursive terms. We derive and analyze the solution of this equation for the case that the αi(z)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _i(z)$$\end{document} are commutative contraction mappings. The results are applied to a wide range of queueing, autoregressive and branching processes.


Introduction
In many applied probability models, in particular in queueing models, the following type of recursion describes the behavior of a key performance measure: where the involved random variables are nonnegative and integer valued. Under certain independence assumptions and ergodicity conditions, this gives rise to the following type of functional equation for the probability generating function (pgf) f (z) of the steady-state distribution of the X n process: An example of such a model is a branching process with immigration (BPI). In that case, the function f (z) represents the pgf of the steady-state distribution of the number of individuals and the function α(z) the pgf of the offspring distribution. If in all states except state 0 the pgf of the immigration distribution is given by g(z) while in the special state 0 the pgf of the immigration distribution is given by g 0 (z), then the function f (z) satisfies (1) with K (z) = f (0) (g 0 (z) − g(z)). Remark that in the special case that g 0 (z) = g(z), we have that K (z) = 0. The solution of (1) is given by an expression containing an infinite product and an infinite sum (see Eq. (13) later on in this paper) which is obtained after iteration of Eq. (1). Branching processes with immigration appear for example in the analysis of the M/G/1 queue with (single or multiple) gated vacations (see Takagi [18], Sect. 2.5 of Chapter 2), the M/G/1 queue with permanent customers (see Boxma and Cohen [3]) and, a multi-type variant, in the analysis of polling systems (see Resing [15]).
In the case that the branching process with immigration evolves in an i.i.d. random environment (BPIRE) in which the environment can be in M different states, the corresponding functional equation is of the form where p i is the probability that the environment is in state i and α i (z) is the pgf of the offspring distribution when the environment is in state i, i = 1, 2, . . . , M.
In this paper, we obtain the solution of functional equation (2) in the particular case that the functions α 1 (z), . . . , α M (z) are commutative contraction mappings on the closed unit disk. Although (2) with multiple recursive terms is a natural extension of (1) with only a single recursive term, it is hardly studied in the queueing literature, probably because the number of different terms after the nth iteration of (2) grows exponentially. That is also the reason why in this paper we restrict ourselves to the case in which the functions α 1 (z), . . . , α M (z) are commutative.
In Adan et al. [1], a specific example of (2) was analyzed in detail in the study of a queueing system with two classes of impatient customers. The main goal of the present paper is to give a general treatment of Eq. (2) and to show how in the commutative case the growth of the number of iteration terms can be handled. An additional aim is to treat several queueing and branching-type examples in which (2) appears, also allowing complications like the occurrence of a pole in g(z) and K (z).
Organization of the paper In Sect. 2, we solve the recursion (2), both for the homogeneous case where K (z) ≡ 0 (Sect. 2.1) and for the inhomogeneous case (Sect. 2.2). The results are applied in the subsequent sections. We start, in Sect. 3, with a partic-ular queueing model. Its choice is motivated by the fact that it provides a relatively simple illustration of the theory for the homogeneous case, while still having a few features that deviate from the setting of Sect. 2. Section 4 considers a special case of the branching process with immigration in a random environment, also called random coefficient integer-valued autoregressive process of order 1, in which the offspring of an individual in each environmental state can only be equal to 0 or 1. In Sect. 5, we consider an integer-valued reflected autoregressive process, which may be viewed as a generalization of an embedded queue length process in the M/G/1 queue. Section 6 is devoted to another reflected autoregressive process, this time on [0, ∞). Some topics for further research are mentioned in Sect. 7.

The recursion
In this section, we study recursion (2) for the generating function f (z) of a nonnegative discrete random variable X with E[X ] < ∞, where g(z) and K (z) are analytic functions (and not necessarily generating functions), p 1 , . . . , p M is a probability distribution and α 1 (z), . . . , α M (z) are commutative contraction mappings on the closed unit disk, i.e., there is a constant κ < 1 such that |α i (z) − α i (u)| ≤ κ|z − u| and α i (α j (z)) = α j (α i (z)) for each i and j. For example, the mappings α i (z) = 1−a i +a i z with |a i | < 1 are contractions with κ = max(a 1 , . . . , a M ), and they commute, since Note that the commutativity property implies that the contractions α i (z) have the same fixed point a. In the example above, we have a = 1. Equation (2) is suitable to iteratively determine f (z). We distinguish between the following two cases.

The homogeneous case K(z) = 0
After n iterations of the homogeneous equation we obtain where α i 1 ,...,i M (z) = α i 1 1 (α i 2 2 (· · · (α i M M (z)) · · · )) and α n i (z) is defined as the nth iterate of α i (z). In particular, α 0,...,0 (z) = z. The functions L i 1 ,...,i M (z) are recursively defined by  (3,2). Its path weight is the product of the step weights g(α 0,0 (z))g(α 1,0 (z))g(α 1,1 (z))g(α 2,1 (z))g(α 3,1 (z))  (4), we proceed as follows. Note that, when i 1 + · · · + i M = n, and hence, for |z| ≤ 1, where the integral is along the segment connecting α i 1 ,...,i M (z) with a. Next we rewrite (4) as follows For given z, the second term in the right-hand side of (8) converges to zero for n → ∞. This can be seen as follows. By substituting z = a in (3), we get g(a) = 1. Hence, from (6), step weight g( In other words, the weight of sufficiently long paths grows at most with 1 + per step. So there is a constant C such that for all n, the weight of a path from (0, . . Then, from (7), we can conclude that the second term in (8) is , which goes to zero for n → ∞. We have thus proven the following theorem.

Theorem 1 (Homogeneous case) The generating function f (z) is given by
Remark 1 The unknown f (a) follows by substituting z = 1 in (9) and using f (1) = 1. This gives

Remark 2
The first term in the right-hand side of (8) is a sum of M+n M−1 terms. This number is O(n M−1 ), which grows quickly (polynomially) in n for already moderate values of M. However, this is not as quickly as in the non-commutative case, in which case we would have M n+1 terms, growing exponentially in n.

Remark 3
Above we have seen that the second term in the right-hand side of (8) converges to zero geometrically fast (with rate (1 + )κ). Hence, the first term in the right-hand side of (8) will provide an accurate approximation for f (z) already for small values of n.

The inhomogeneous case K(z) = 0
We now consider inhomogeneous Eq. (2) and obtain after n iterations, The first term of (10) is the same as in the homogeneous case. For given z, we proved that it converges if g(a) = 1. For convergence of the second term, we need that either In this subsection, we successively consider both cases.
, in this case we basically assume that K (a) = 0. Note that, by substituting z = a in (2), this assumption implies that g(a) = 1 if f (a) = 0 and thus that the first term in (10) converges. We now show that it also implies convergence of the second term. Similar to (7), we have for i 1 + · · · + i M = n, where D is the maximum value of |K (u)| in the closed disk with center a and radius |z − a| (intersected with the closed unit disk). Then, the kth term in the double sum appearing in (10) is bounded by C(1 + ) k κ k |z − a|D, and hence, the double sum converges for n → ∞. We conclude that the following holds.
Theorem 2 (Inhomogeneous case) Provided K (a) = 0 and f (a) = 0, the generating function f (z) is given by
Theorem 3 (Inhomogeneous case) Provided |g(a)| < 1, the generating function f (z) is given by Remark 5 In this section, we studied Eq. (2) for a generating function f (z) of a nonnegative discrete random variable X . It is readily seen that Theorems 1-3 are also valid in case f (z) is the Laplace-Stieltjes transform (LST) of a non-negative continuous random variable X . Of course, then α 1 (z), . . . , α M (z) are assumed to be contraction mappings on the closed positive half space (instead of the closed unit disk).

The D M /G/1 shot-noise queue
In this section, we consider the workload at arrival epochs for a specific queue. The analysis of the LST of the workload gives rise to a simple recursion that has the form (2), in the homogeneous variant. This section thus provides a relatively simple illustration of our theory for the homogeneous case-while still having a few features that deviate from the setting of Sect. 2. The model under consideration is the D M /G/1 shot-noise queue; we refer to [7] for a recent survey on shot-noise queueing models. The D M /G/1 shot-noise queue is a single server queue, in which the successive interarrival times A 1 , A 2 , . . . of customers are i.i.d., with the distribution of a generic interarrival time A being given by The service requirements of successive customers B 1 , B 2 , . . . are i.i.d. with finite mean, and with LST β(·); all interarrival times and service times are independent. The special feature of the model is that the server speed is workload proportional (shot noise): when the workload is x, the service speed is r x. Let X n denote the workload just before the arrival of the nth customer. It is well known [2] that, in between arrivals, the workload decreases exponentially; hence, It is readily verified that, due to the workload-proportional decrease, the steady-state distribution of the {X n , n = 1, 2, . . . } process exists. Stability conditions for queueing and storage models with more general workload-dependent decay have been discussed, a.o., by Brockwell et al. [10] and Cinlar and Pinsky [11] . Let X denote a random variable distributed as the steady-state distribution of the process {X n , n = 0, 1, . . . }, with LST ξ(·), then ξ(s) = β(as)ξ(as) = β(as)β(a 2 s)ξ(a 2 s) = · · · = ξ(a n s) Now observe the following. Firstly, ξ(a n s) → ξ(0) = 1 for n → ∞. Secondly, convergence is geometric, as |1 − ξ(a n s)| ≤ ∞ 0 |1 − e −a n st |dP(X < t) ≤ a n s ∞ 0 tdP(X < t) = a n sE [X ]. [19]); hence, the convergence of the product in (18) follows from We conclude that Some thought will make it clear that the ith term in this infinite product represents the contribution to X from an arrival that occurred i arrivals before the present one.
Let us now turn to the general D M /G/1 shot-noise case, cf. (15). After n iterations, (17) gives (very similarly to the analysis in Sect. 2.1) Notice that an (i 1 , . . . , i M ) term corresponds to a contribution to the workload (just before an arrival epoch) from an arrival that occurred i 1 + · · · + i M arrivals before the present one, the total interval consisting of i k interarrival times of length t k , k = 1, . . . , M, in any of i 1 +···+i M i 1 ,...,i M orders. It is readily seen that and hence, the sum in (22)  we have the following theorem.

Theorem 4
The LST of the steady-state workload just before arrival epochs in the D M /G/1 shot-noise queue is given by

Remark 6
One could subsequently derive the steady-state workload LST at an arbitrary epoch by averaging over one arrival interval, and using a stochastic mean-value theorem.
where the Z n are nonnegative integer-valued random variables with finite mean. The Y k,n are Bernoulli random variables with parameter ξ n , i.e., P(Y k,n = 0) = 1 − ξ n and P(Y k,n = 1) = ξ n ; but we assume the special feature that the ξ n are themselves random variables, independent and identically distributed with P(ξ n = a i ) = p i , i = 1, . . . , M. Hence, in generation n it holds with probability p i , i = 1, . . . , M, for all Y k,n that they are 1 with probability a i and 0 with probability 1 − a i . All Z j and Y k,m are also assumed to be independent. In Sect. 4.1, we consider the steadystate distribution of the process {X n , n = 0, 1, . . . }, and in Sect. 4.2 we do this for the generalization in which the BPIRE process behaves differently at zero, i.e., when X n = 0 for some n.

The steady-state case
The generating function, f (z), of the stationary distribution of the process {X n , n = 0, 1, . . . } satisfies the recursion where g(z) is the pgf of the random variable Z n . Hence, we are in the homogeneous case (3) with contraction mappings α i (z) = 1 − a i + a i z. In this case, the functions α i 1 ,...,i M (z) = α i 1 1 (α i 2 2 (· · · (α i M M (z)) · · · )) are given by Define the functions L i 1 ,...,i M (z) again recursively by (5) and use that the contraction mappings α i (z) have fixed point a = 1 in this case, and hence, f (a) = f (1) = 1. Theorem 1 now implies the following theorem.

Theorem 5
The steady-state probability generating function f (z) of the BPIRE process is given by

Remark 7
In the special case that p 1 = 1, the recursion becomes yielding In the special case that a 2 = · · · = a M = 0, the recursion becomes and in this case the solution is given by

Deviating behavior at zero
In this subsection, we assume that the BPIRE process behaves differently at zero, i.e., when X n = 0 for some n. In particular, we assume that with V 0 , V 1 , . . . i.i.d. nonnegative integer random variables, with pgf g 0 (z). V n is assumed to be independent of all Y i,m , Z m and X m , m = 0, 1, . . . , n. It is readily seen that the steady-state pgf f (z) in this case satisfies the recursion Hence, we are in the inhomogeneous case of Eq.
For future use we observe, by substituting z = 0 in (37), that We conclude that the following holds.

A reflected autoregressive process and an M/G/1 queue generalization
In this section, we consider an integer-valued stochastic process {X n , n = 0, 1, . . . } that is very similar to the autoregressive process of the previous section, but includes a negative component and has reflection at zero. It is determined by the following relation: with [x] + = max(0, x), Z 1 , Z 2 , . . . i.i.d. nonnegative integer-valued random variables with pgf C(z) and where Y k,n are i.i.d. Bernoulli distributed random variables: P(Y k,n = 1) = ξ n , P(Y k,n = 0) = 1 − ξ n . In addition, we (again) assume the special feature that the ξ n are themselves random variables, independent and identically distributed with Hence, in generation n it holds with probability p i , i = 1, . . . , M, for all Y k,n that they are 1 with probability a i and 0 with probability 1 − a i . All Z j and Y k,m are also assumed to be independent of each other and of all preceding X r .
If all Y k,n are equal to one, then (40) can be interpreted as follows. Consider the number of waiting customers in the M/G/1 queue, just after the beginning of the nth service. Let X n denote this number, and let Z n denote the number of arrivals during the nth service. Then, X n+1 = [X n + Z n − 1] + . If, in addition, each of the X n customers becomes impatient with probability 1 − a during the nth service and leaves, then the sequence {X n } satisfies (40) with p 1 = 1 and a 1 = a, i.e., with ξ n ≡ a.
Below we show how the pgf f (z) of the steady-state distribution of the {X n , n = 0, 1, . . . } process can be obtained. It follows from (40), with [x] − = min(0, x) and X denoting a generic random variable with pgf f (z), that (44) In [8], both the transient behavior and steady-state behavior of the R n process with A n ≡ a are studied, via a Wiener-Hopf technique (cf. [12]) that leads to a recursion. We apply the same tools in the extension defined by (43), (44). Below we first follow the approach of [8]. Introduce U n := min(A n R n + G n , 0) for n = 0, 1, . . . , and the transforms The first transform is analytic for Re s ≥ 0 and the second one for Re s ≤ 0. Observing that 1 + e x = e max(x,0) + e min(x,0) , we have for n = 0, 1, . . . : e −s R n+1 = e −s(A n R n +G n ) + 1 − e −sU n .
Taking expectations and realizing that R n , A n and G n are independent, we can write E e −s R n+1 |R 0 = z = E e −sG n M i=1 p i E e −sa i R n |R 0 = z + 1 − E e −sU n |R 0 = z , (46) hence, after multiplication of both sides by r n+1 and summation, we obtain for Re s = 0: Restricting ourselves at this stage to Re s = 0 ensures that all terms are properly defined. Multiplying both sides by λ − s, one obtains: Because all a i < 1, the steady-state distribution of the {R n , n = 0, 1, . . . } process always exists [13]. We shall restrict ourselves to the steady-state case. (The transient case can in principle be studied in a similar way. Here it should be observed that, with fixed point a = 0, we have K (0) = 1 + r 1−r − rU z (r , 0) = 1 = 0. We are now in Case 2 of Sect. 2.2; |r | < 1 will guarantee the convergence of the corresponding p i 1 1 . . . p i M M L i 1 ,...,i M (z) as i 1 + · · · + i M = n → ∞.)