Some reflected autoregressive processes with dependencies

Motivated by queueing applications, we study various reflected autoregressive processes with dependencies. Amongst others, we study cases where the interarrival and service times are proportionally dependent with additive and/or subtracting delay, as well as cases cases where interarrival times depends on whether the service duration of the previous arrival exceeds or not a random threshold. Moreover, we study cases where the autoregressive parameter is constant as well as a discrete or a continuous random variable, as well as cases where . More general dependence structures are also discussed. Our primary aim is to investigate a broad class of recursions of autoregressive type for which several independence assumptions are lifted, and for which a detailed exact analysis is provided. We provide expressions for the Laplace transform of the waiting time of a customer in the system in terms of an infinite product of known Laplace transforms. An integer-valued reflected autoregressive process that can be used to model a novel retrial queueing system with orbital searching time to depend on whether the last busy period starts with an empty or a non empty orbit queue, is also discussed. For such a model the probability generating function of the stationary orbit queue length is given as an infinite product of known generating functions. A first attempt towards multidimensional setting is also analyzed. Some additional generalizations with more general dependence structure are also discussed.


Introduction
This work focuses on various stochastic recursions of autoregressive type, such as: , with probability (w.p.) q := 1 − p.
The ultimate goal of this work is to investigate classes of reflected autoregressive processes described by recursions of the type given above, in which various independence assumptions of {B n } n∈N0 , {A n } n∈N0 are lifted and for which, a detailed exact analysis can be also provided.
The stochastic recursion (1) where {V n } n∈N0 are constant and where {B n } n∈N0 , {A n } n∈N0 are i.i.d.sequences, and also independent on {W n } n∈N0 has been treated in [7], i.e., the case where W n+1 = [aW n + B n − A n ] + , n = 0, 1, . .., with a ∈ (0, 1).The case where a = 1 corresponds to the classical Lindley recursion describing the waiting time of the classical G/G/1 queue, while the case where a = −1 is covered in [12].Further progress has been made in [5], where additional models described by recursion (1) have been investigated.The work in [5,Section 3] is the close to our case, where the authors investigated a recursion where V is either a positive constant with probability p, or a random variable taking negative values with probability 1 − p.The fact that V is negative simplified considerably the analysis.
In [4], the authors considered the case where V n W n in (1) was replaced by S(W n ), where {S(t)} is a Levy subordinator (recovering also the case in [7], where S(t) = at).Note that in [7,5,4] the sequences {B n } n∈N0 , {A n } n∈N0 are assumed to be independent.Recently, in [6], the authors considered Lindley type recursions that arise in queueing and insurance risk models, where the sequences {B n } n∈N0 , {A n } n∈N0 obey a semi-linear dependence.These recursions can also be treated as of autoregressive type.This work is the closer to ours.Moreover, in [1], the authors developed a method to study functional equations that arise in a wide range of queueing, autoregressive and branching processes.Finally, the author in [9], considered a generalized version of the model in [5], by assuming V n to take values in (−∞, 1].In particular, in [9], the author investigated the recursion (4) for a = 1.
The main contribution of this paper is to investigate the transient as well as the stationary behaviour of a wide range of autoregressive processes described in (1)- (5), by lifting various independence assumptions of the sequences {B n } n∈N0 , {A n } n∈N0 .This is accomplished by using Liouville's theorem, and by stating and solving a Wiener-Hopf boundary value problem, or by solving an integral equation, depending on the nature of {V n } n∈N0 .We have to point out here that to our best knowledge, autoregressive recursions of the form (2)- (5) have not considered in the literature so far.We also investigate the stationary analysis of {X n } n∈N0 in (5), which represents a novel retrial queueing model.An extension to a two-dimensional case that describes the priority case is also considered.

M/G/1-type autoregressive queues with interarrival times randomly proportional to service/system times
In the following we cope with some autoregressive M/G/1-type queueing systems where the interarrival time between the nth and the (n + 1)th job, say A n , depends on the service time of the nth job, or on the system time after the arrival of the nth job.

Interarrival times randomly proportional to service times
Consider the following variant of the standard autoregressive M/G/1 queue: When the service time equals x ≥ 0 following the nth arrival, then the next interarrival time equals β i x (with probability p i , i = 1, . . ., N + M ) increased by an independent additive delay J n .In the following, we consider the recursion (1), where P (V n = a) = 1, a ∈ (0, 1).Let W n the waiting time of the nth customer.The interarrival time between the nth and the (n + 1)th customer, say A n , satisfies A n = G n B n + J n , where B n is the service time of the nth customer and J n an additive delay or random jitter.The random variable G n has finite support.Let β i denote its ith value and let p i = P (G n = β i ) denote the corresponding probability, i = 1, . . ., N + M (M, N ≥ 1), with N +M i=1 p i = 1.We further assume that the service times and jitter are exponentially distributed: B n ∼ exp(µ) and J n ∼ exp(δ).Extensions to the case where J n has a rational transform will be also discussed.Thus, the sequence {W n } n∈N0 obeys the following recursion: where a ∈ (0, 1).Without loss of generality and in order to avoid trivial solutions, assume that 1 < β 1 < β 2 < . . .< β N , and β N +1 < β N +1 < . . .< β N +M < 1.

Transient analysis
We first focus on the transient distribution, and following the lines in [7], let for |r| < 1, where Then, using the property that 1 + e x = e [x] + + e Assume hereon that B ∼ exp(µ).Then, φ B ( βi s) = µ µ+ βis = 1 1−γis , where Note that g(s) = 0 has N + M distinct and real roots γ −1 i , where N of them are positive and M are negative.In particular, let s . Now (7) becomes for Re(s) = 0: Now we make the following observations: • The left-hand side is analytic in Re(s) > 0, and continuous in Re(s) ≥ 0, • The right-hand side is analytic in Re(s) < 0, and continuous in Re(s) ≤ 0, Thus, (8) represents an entire function.Generalized Liouville's theorem [11,Theorem 10.52] states that in their respective half planes, both the left and the right-hand side can be written as the same (N + 1)−th order polynomial in s, for large s, i.e., Note that for s = 0 (9) yields Having in mind that f (0) , we easily realize that C 0.w (r) = 0.Moreover, setting s = δ, and s = s + i , i = 1, . . ., N we obtain the following system of equations for the rest of unknowns C i,w (r), i = 1, . . ., N : −rδ It remains to obtain Z w (r, as + i ), i = 1, . . ., N , Z w (r, aδ).These terms are evaluated as follows: Expression ( 9) is now written as Z w (r, s) = K w (r, s)Z w (r, as) + L w (r, s), where Iterating (11) yields with the convention of empty product to be equal to one.Setting s = αδ, and s = as + i in , we obtain expressions for the Z w (r, aδ), Z w (r, as + i ), i = 1, . . ., N , respectively.Substituting back in (10), we obtain an (N +1)×(N +1) system of equations with C i,w (r), i = 1, . . ., N + 1 as unknowns.
Remark 1 It is easily realised in (12) that Z w (r, s) appears to have singularities in s = δ/a j , and s = s + i /a j , i = 1, . . ., N , j = 0, 1, . ... We can show that these are removable singularities.Let us show this for s = δ, and s = s + i .We write (12) as follows to isolate the singularities for s = δ, and s = s + i : It is easily realised by using (10), that the term inside the brackets in the last line vanishes for s = δ, and s = s + i , confirming that s = δ, and s = s + i , i = 1, . . ., N are not poles of Z w (r, s).Similarly, we can show using (9) that Z w (r, s) has no singularity at s = δ/a, s = s + i /a, and so on.

Stationary analysis I
One can follow an analogous analysis as the one presented in the previous sub-subsection or the lines in [7] and obtain a functional equation similar to (11).In particular, by applying Abel's theorem, or by considering the relation Since a ∈ (0, 1), the stability condition can be ensured as long as E(log(1 + (1 − G)B)) < ∞ (see also [5]).Note also that where we used the fact that E(e −s(aW∞+ is the LST of the pdf characterized by: and thus, E(e −s(aW∞+ 13) is now written as since lim n→∞ Z(a n s) = Z(0) = 1.Note that P = P (W ∞ = 0), and is derived by multiplying (14) with δ − s (i.e., the functional equation before the iterations), and setting s = δ, so that .

Stationary analysis II
We now focus on the stationary distribution of W n as n → ∞: In the following, we follow a different approach compared to the previous sub-subsection: Denote by Z(s) := lim n→∞ Z n (s), where Z n (s) := E(e −sWn ) and let f Wn (w) the probability density function (pdf) of W n (for convenience, assume also W 0 = 0).Assume now that β i ∈ (0, 1), i = 1, . . ., K, where K = N + M (so that Re(s βi ) ≥ 0 for Re(s) ≥ 0), and service time is arbitrarily distributed with pdf f B (.), and LST φ B (.).Then, (15) Letting n → ∞ we get: Iterating ( 16) we get: Differentiating ( 16) with respect to s, and setting s = 0 yields after some algebra where Z(aδ) as given in (17), and E(B) the mean service duration.
Remark 2 Note that the analysis can be considerably adapted to consider the case where J n follows a hyperexponential distribution with L phases, i.e., with density function f J (x) : where where .

Proportional dependency with additive and subtracting delay
We now focus on the case where the interarrival times are such that A n = cB n + J n , with J n := Jn , with probability p, − J n , with probability q := 1 − p, where Jn ∼ exp(δ), J n ∼ exp(ν).Now the sequence ( With probability p, J n = Jn , and thus (cB n + J n ) + = cB n + Jn , while with probability q, J n = J n , and thus where c := 1 − c.Assume that the steady-state waiting time distribution exists and denote by Z(s) its LST.Following similar steps as in the derivation of (15), we can obtain: H(a j s).
Setting s = aδ, and substitute back, we obtain Remark 4 One may also consider the case where the interarrival times are related to the previous service time as follows: A n = G n B n + J n , where J n , as given in (22), and G n are iid random variables with probability mass function given by and following the same arguments as above, we again have where now Following the lines in subsection 2.1.2,and having in mind that lim n→∞ Z(a n s) = 1, we obtain the desired expression for Z(s).

Remark 5
The case where Jn , J n are i.i.d.random variables following a distribution with rational LST can also be treated similarly.In particular, assume that Jn , J n follow hyperexponential distributions, i.e., their pdfs are f J (x) := L i=1 q i δ i e −δix , and f J (x) := M i=1 h i ν i e −νix , respectively.Then, following similar arguments as above, and assuming A n = G n B n + J n , where J n , as given in (22), we obtain after lengthy computations: where now Iterating (25) as in subsection 2.1.2,and having in mind that lim n→∞ Z(a n s) = 1, we obtain the desired expression for Z(s).
Note that the recursion (26) is a special case of the recursion (1) with V n := 1 − G n .

Asymptotic expansions
In the following we focus on deriving asymptotic expansions of the basic performance metrics P (W = 0), E(W l ), l = 1, 2, . .., by perturbating β i s, i.e., by letting in (27) β i to be equal to β i ǫ with ǫ very small.The, ( 27) is written: Note that for ǫ = 0, the above equation provides the LST of the waiting time (say W ) of the classical M/G/1 queue where arrivals occur according to a Poisson process with rate δ.So when ǫ → 0, there is a weak dependence between soji=ourn time and the subsequent interarrival time.Following [6, subsection 2.3], consider the Taylor series development of P (W = 0), E(W l ), l = 1, . . ., L up to ǫ m terms for m ∈ N.So for ǫ → 0: Differentiating the functional equation with respect to s, setting s = 0 yields for ρ = δE(B) Assuming that the first L moments of W are well defined, we subsequently differentiate the functional equation above l = 2, . . ., L times with respect to s, and set s = 0.Then, for l = 2, 3, . . ., L we have: Setting ǫ = 0, and having in mind that K i=1 p i = 1 we recover the recursive formula to obtain the moments of the standard M/G/1 queue: Then, substituting (28) in (29) we have (30) Equating ǫ factors on both sides we obtain R l−1,1 in terms of R l−2,1 , . . ., R 0,1 , as well as in terms of E( W n ) obtained above.Thus, since R 0,1 is known, all R l,1 can be derived by: Similarly, for h = 2, The procedure we follow to recursively obtain R l,h is the same as the one given in [6, subsection 2.3], so further details are omitted.

The single-server queue with service time randomly dependent on waiting time
Customers arrive according to a Poisson process with rate λ.If the waiting time W n of the nth arriving customer, then her service time equals [B n −Ω n W n ] + , with P (Ω n = a l ) = g l , l = 1, . . ., K, with a l ∈ (0, 1), where {B n } n∈N0 is a sequence of independent, exponentially distributed random variables with rate µ, independent of anything else.Note that when the waiting time is very large the service requirement tends to zero, which can explained as an abandonment.Let Z(s) := E(e −sWn ), and A n i.i.d.random variables from exp(λ).Then, where Assuming n → ∞ and studying the limiting variable W we have Thus, substituting back in (31) we arrive after simple calculations in: where . Note also that Z(µa l ) = P (B > a l W ), and K l=1 g l Z(µa l ) = P (B > ΩW ).We need to iterate (32) and having in mind that as s → ∞, Z(s) → 0 (needs some work).Note that such kind of recursions were treated in [1], since the cummutativity of ζ l (s) := s + a l µ and ζ m (s ) makes the recursion (32) relatively easy to handle, although in each iteration, any term gives rise to new K terms; see also [6,Remark 5.3].Extensions to the case where service time distributions have rational LST are doable; e.g.hyperexponential.
4 Threshold-type dependence among interarrival and service times

The simple case
Customers arrive with a service request at a single server.Service requests of successive customers are independent, identically distributed (i.i.d.) random variables B i , i = 1, 2, . . .with distribution B(.), mean b and Laplace-Stieltjes transform (LST ) β(.).Upon arrival, the service request is registered.If the service request B i is less than a threshold T i , then the next interarrival interval, say J (0) i , is exponentially distributed with rate λ 0 ; otherwise, the service time becomes exactly equal to T i (is cut off at T i ), and the next interarrival interval, say i , is exponentially distributed with rate λ 1 .We assume that an arrival makes obsolete a fixed fraction 1 − a 0 (resp. 1 − a 1 ) of the work that is already present, with a k ∈ (0, 1), k = 0, 1.We assume the threshold T i to be i.i.d.random variables with general distribution T (•), that has mean τ and LST τ (•).Let also for Re(s) ≥ 0 Let W n the waiting time of the nth arriving customer, n = 1, 2, . ...Then, where Then, (34) leads for Re(s) = 0 to: Multiplying (35) by Our objective is to obtain Z w (r, s), M w (r, s) by formulating and solving a Wiener-Hopf boundary value problem.
A few observations: • The LHS in (36) is analytic in Re(s) > 0, and continuous in Re(s) ≥ 0.
• Z w (r, s) is for Re(s) ≥ 0 bounded by | r 1−r |, so by the generalized Liouville's theorem [11,Theorem 10.52], the LHS is at most a quadratic polynomial in s (dependent on r) for large s, Re(s) > 0.
• M w (r, s) is for Re(s) ≤ 0 bounded by | r 1−r |, so by the generalized Liouville's theorem [11,Theorem 10.52], the RHS is at most a quadratic polynomial in s (dependent on r) for large s, Re(s) < 0. Thus, with C i,w (r), i = 0, 1, 2, functions of r to be determined.Taking s = 0 in (37) yields and having in mind that χ(0 Thus, C 1,w (r), C 2,w (r), needs to be determined, since Z w (r, α i λ i ), i = 0, 1 are still unknown.This difficulty can be overcome by deriving an expression for Z w (r, α i λ i ), i = 0, 1, by deriving Z w (r, s) after successive iterations of (37), which now becomes where, After n − 1 iterations, we obtain where K k,n−k (s) are recursively defined as follows: Therefore, The second term in the RHS of (44) converges to zero due to the fact that |r| < 1, thus, Setting in (45) s = a i λ i we obtain expressions for the Z w (r, a i λ i ), i = 0, 1.Note that these expressions are given in terms of the unknowns C i,w (r), i = 0, 1. Substituting back in (39), (40), we obtain a linear system of two equations with two unknowns C i,w (r), i = 0, 1.
Differentiating (50) with respect to s and letting s = 0 yields after some algebra that, where f ′ (.) denotes the derivative of a function f (.) and C 1 , C 2 are derived as shown above.

4.1.2
The case a 0 ∈ (0, 1), a 1 = 1 We now consider the stationary version of the special case where a 1 = 1, i.e., we assume that when B n ≥ T n , the next arrival does not make obsolete a fixed fraction of the already present work.This can be considered natural if we think that in such a case the service time is cut-off since exceeds the threshold T n .Following similar arguments as above, we obtain where β(s) := χ(s) is the LT of the distribution of the random variable B, which is the time elapsed from the epoch a service request arrives until the epoch the registered service is of threshold type: .
Thus, following the lines in [7], Liouville's theorem [11] states that For s = 0, (55) implies that C 0 = 0. Thus, which has a solution similar to the one in [7, Theorem 2.2], so further details are omitted.

4.1.3
The case a 0 = a 1 := a ∈ (0, 1) Now consider the case where the fraction of work that becomes obsolete because of an arrival, is independent on whether B < T or B ≥ T .In such a scenario, • The LHS of (57) is analytic in Re(s) > 0, and continuous in Re(s) ≥ 0.
• Z(s) is for Re(s) ≥ 0 bounded by 1, and hence the LHS of (57) behaves at most as a quadratic polynomial in s fro large s, with Re(s) > 0.
• M (s) is for Re(s) ≤ 0 bounded by 1, and hence the RHS of (57) behaves at most as a quadratic polynomial in s fro large s, with Re(s) < 0.
Remark 9 Assume now that the interarrival times deterministic proportionally dependent to service times.More precisely, let n := T n , and Following similar arguments as in the previous section, we arrive for Re(s) = 0, at Using similar arguments as above, Liouville's theorem [11] implies that The rest of the analysis follows as the one in the previous section.Similar steps as those in the previous section can be followed to cope with the stationary analysis, so further details are omitted.
• For large s, both sides in (69) are O(s L0+L1 ) in their respective half-planes.

A more general case
We now consider the case where the interarrival times are also depended on the system time.More precisely, we assume that J Then, by using (75) and using similar arguments as above we obtain for Re(s) = 0, Z w (r, s) − re −sw = rφ X0 (−s) K i=1 p i χ(s βi )Z w (r, s βi ) + rφ X1 (−s) M i=1 q i ψ(sγ i )Z w (r, sγ i ) + r 2 1−r − rM w (r, s), or equivalently, Now, • The LHS of (76) is analytic in Re(s) > 0, and continuous in Re(s) ≥ 0.
• For large s, both sides in (69) are O(s L0+L1 ) in their respective half-planes.

A mixed case
Consider the following recursion: where V n < 0, and a ∈ (0, 1).Then, following a similar procedure we obtain Equivalently, we have Clearly, • the LHS of (84) is analytic in Re(s) > 0 and continuous in Re(s) ≥ 0, • the RHS of (84) is analytic in Re(s) < 0 and continuous in Re(s) ≤ 0, • for large s, both sides are O(s 2 ) in their respective half planes.

The uniform proportional case with dependence
In the following, we consider recursions of the form with V n ∼ U (0, 1), and dependence among B n , A n .
Let U W (r, s) := ∞ n=0 r n Z n (s), |r| < 1, then: The solution of such kind of first-order differential equation is obtained by following the lines in [5, Section 5], so further details are omitted.

Randomly proportional dependency with additive delay
In the following we consider the case where . ., K, and J n to follow a hyperexponential distribution with pdf f (x) = L j=1 q j δ j e −δj x (the analysis can be further generalized in the case of a distribution with a rational Laplace transform).Then, e −s(vw+ βl x−y) δ j e −δjy + ∞ y=vw+ βl x δ j e −δj y dy dxdP (W < w)dv where Q the type of the arrival process.Then, multiplying with r n+1 and sum over n (with W 0 = w) results in where W (r, s), j = 1, . . ., L. Then, where The form of ( 94) is that same as in [4, Section 3], and the analysis can be performed similarly.

Proportionally dependent on system time
We now consider the case where A n = c(W n + B n ) + J n , c ∈ (0, 1).We assume that (B n , J n ) are i.i.d.sequences of random vectors.Thus, the quantities (cB n − J n ) are i.i.d.random variables, however, within a pair B n , J n are dependent.Here, we assume that a non-negative random vector (B, J) has a bivariate matrixexponential distribution with LST E(e −sB−zJ ) := G(s,z) D(s,z) , where G(s, z), D(s, z) are polynomial functions in s, z.A consequence of this definition is that the transform of Y := cB − J is also a rational function; the distribution of Y is called a bilateral matrix-exponential [2,Theorem 3.1].This class of distributions, under which we model the dependence structure, belongs to the class of multivariate matrix-exponential distributions, which was introduced in [3].For ease of notation, let E(e −sY ) := h(s) = f (s) g(s) , and assume that g(s) has L zeros, say t j such that Re(t j ) > 0, j = 1, . . ., L, and M zeros, say h j , such that Re(h j ) < 0, j = 1, . . ., M , whereas f (s) is a polynomial of degree at most M + L − 1, not sharing the same zeros with g(s).
Note that now recursion (92) becomes Multiplying with r n+1 (0 < r < 1) and summing for n = 0 to infinity, and having in mind that W 0 = w, we obtain where H(r, s) = ∞ n=0 r n E(e −sHn ).Now we have: 1. the LHS of (95) is analytic in Re(s) > 0 and continuous in Re(s) ≥ 0, 2. the RHS of (95) is analytic in Re(s) < 0 and continuous in Re(s) ≤ 0, 3. for large s, both sides are O(s M+L ) in their respective half planes.
Since we will let n to tend to ∞ we are interested in investigating the convergence of the summation in the previous expression as well as in obtaining the limit of the first term in the right hand-side of the previous expression.Since the expressions of K(r, s), L(r, s) share the same properties as those in [9], we can show that We still need M more equations to obtain a system for the coefficients C l (r).Substituting s = h j , j = 1, . . ., L in (97) and using (100) we obtain Finally, by using ( 98), (101) we can derive the remaining coefficients C l (r), l = 1. . . ., L + M .
Remark 10 An alternative way to solve (99) is by performing the transformation v 1 = sy 1 , so that (99) becomes: Note that (102) is a Fredholm equation [10]; therefore, a natural way to proceed is by successive substitutions.Define now iteratively the function with L 0 * (r, s) := L(r, s).Then, after n iterations we have that Note that To see this, observe that Thus, the above limit is less than or equal to Note that from (106) w (r, u)du, Assume now that A Thus, following standard techniques from the theory of ordinary differential equations, we have for a positive number c, such that c < s, Since I (a)′ (s) = s a−1 Zw(r, s), I (a)′ (0) = 0, and thus, Combining the above with (109), and having in mind that Now substituting the derived expression for Z w (r, s) in (107) we can obtain a system of equations to obtain the remaining unknowns c k (r), k = 1, . . ., l.

Other generalizations
We consider the case where , and E(e −sAn |B n = t) = χ(s)e −ψ(s)t , and B n ∼ exp(µ).Thus, the interarrival times distribution depends on the service time of the previous customer, so that when Re(µ + ψ(s) + z) > 0. Since for s = 0 the E(e −sAn |B n = t) should be equal to one, we have to implicitly assume that ψ(0) = 0 and χ(0) = 1.Therefore, by denoting Z n (s) = E(e −sWn ) we have: where U − VnWn (s) := E(e −s[VnWn+Bn−An] − ).Clearly, under the transformation v = su, we have: Letting Z w (r, s) = ∞ n=0 r n E(e −sWn ), r ∈ [0, 1) we have for Re(s) = 0 that It is readily seen that • The LHS in (110) is analytic for Re(s) > 0, and continuous for Re(s) ≥ 0, • The RHS in (110) is analytic for Re(s) < 0, and continuous for Re(s) ≤ 0, • For large s, both sides are O(s L+K0+K1+K2 ) in their respective half-planes.
It can easily verified that G k,n−k (r, s) is a sum of n k terms, and each of them is a product of n terms of values of K(r, f .,. ), which are related to the LST φ Y (.).We have to mention that our framework is related to the one developed in [1] with the difference that the functions f i,j (s) (for j > 0) are more complicated compared to the corresponding a i (z) in [1], and inherit difficulties in solving (119).
In what follows we will let n → ∞ in (120) so that to obtain an expression for F (r, s).In doing that, we have to verify the convergence of the summation in the second term in the right hand side of (120), as well as to estimate the limit of the corresponding first term in the right hand side of (120).The key ingredient shall be the boundness of G k,n−k (s).Similarly to [1, p.8], G k,n−k (s) can be interpreted as the total weight of of all n k paths from (0, 0) to (k, n − k).Let C k,n−k the set of all paths leading from (0, 0) to (k, n − k), where a path from (0, 0) to (k, n − k) is defined as a sequence of grid points starting from (0, 0) and ending to (k, n − k) by only taking unit steps (1, 0), (0, 1).Then, a typical term (one of the n k terms) of G k,n−k (s) should be the following: . . .
where Z 1 , Z 2 , . . . is a sequence of i.i.d.random variables with the same distribution as V + .Following the same procedure as in [9, pp. 8-9] we can show that the each of the weights of the path is bounded, implying that G k,n−k (s) is also bounded.This result will imply as n → ∞, that the first term in the right hand side of (120) vanishes, Thus, We are now ready to obtain the coefficients C l (r), l = 1, . . ., M + L. For s = t j , j = 1, . . ., L, in (116), we have (122), (123) serve as a system of equations for the coefficients C l (r), l = 1, . . ., M + L.
Remark 12 Note that in this subsection we have not considered any dependence framework among B n , A n , since our major focus was on introducing this mixed-autoregressive concept, and generalized the work in [9], by assuming a ∈ (0, 1), instead of a = 1.However, the analysis is still applicable even when we lift this assumption.For example, assume the simple scenario where now with probability p 1 , A n = cB n + J n , i.e., the interarrival time is linearly dependent on the service time, with c ∈ (0, 1), J n ∼ exp(δ).Then, (118) becomes Z w (r, sy 1 )µ(dy 1 ) + L w (r, s), where now K 1 (r, s) := r δ δ−s φ B (cs).The rest of analysis can be treated similarly as above.Clearly, the analysis is still applicable if we consider J n to has distribution with rational LST, or more general dependence structure, e.g., A n = G n (W n + B n ) + J n , P (G n = β i ) = q i , i = 1, . . ., M , or the (random) threshold dependence structure analysed in Section 4. Clearly, we can also apply the same steps when lifting independence assumptions for the general case analysed in Remark 11.

A more general dependence framework
In the following we consider a more general dependence structure among {B n } n∈N0 , {A n } n∈N0 .More precisely, assume that thus, the interarrival times distribution depends on the service time of the previous customer, so that with Re(ψ i (s) + z) > 0. Clearly χ(0) = 1, ψ i (0) = 0.The component e −ψi(s)t depends on the previous service time.The component χ(s) does not depend on the service time.
Note that with the above framework we can recover some of the cases analysed above.In particular, the case A n = cB n + J n : N = 1, so that p 1 = 1, E(e −sAn |B n = t) = E(e −s(cBn+Jn) |B n = t) = E(e −sJn )e −cst , with χ(s) := E(e −sJn ), ψ(s) = cs.
The case A n = G n B n + J n , with P (G n = β i ) = p i , i = 1, . . ., N .Then: Another interesting scenario: i=1 H i,k , with probability p i , N i (t) ∼ P oisson(γ i t), and {H i,k } sequences of i.i.d.r.v. with a rational LST, each of them distributed like H i .Then, So, turning back to the simpler general scenario for the stochastic recursion in (1): , where the interarrival times distribution depends on the service time of the previous customer based on (127), we have: s).Assuming that the limit as n → ∞ exists: Assume that χ(s) := , with A 1 (s) a polynomial of order at most K − 1, not sharing the same zeros with the denominator of χ(s), and similarly, B i (s) polynomial of order at most L i − 1, not sharing the same zeros with the denominator of ψ i (s), for i = 1, . . ., N .Then, for Re(s) = 0, By using similar arguments as above, Liouville implies that K j=1 Setting s = 0, yields C 0 = 0.The rest C j s are found by using the K zeros s = λ k , k = 1, . . ., K. Indeed, set s = λ k , k = 1, . . ., K in (128) to obtain the following system: However, we still need to find Z(aλ k ), k = 1, . . ., K.This can be done by iterating This will result in expressions containing infinite products of the form ∞ m=0 A(−a m s)φ B (a m s + ψ i (−a m s)).Indeed, after the iterations we get: (130) Note that for large m, φ B (a m s + ψ i (−a m s)) approaches 1, since a m s + ψ i (−a m s) → 0. Substituting s = aλ k , k = 1, . . ., K in (130), we obtain Z(aλ k ).Finally, by substituting the derived expression in (129), we get a system of equations for the constants C j , j = 1, . . ., K.
Remark 13 Note that in the independent case, i.e., ψ(s) = 0, the situation is easy.In the linear dependent case, i.e., A n = β i B n + J n , ψ i (s) = β i s, the analysis is also easy to handle.If we additionally assume that J n ∼ exp(δ), then, we are interesting in the convergence of ∞ m=0 δφB (a m βis) δ−a m s , which is also easy to handle.
Remark 14 Note that C(z) (resp.G(z)) refers to the pgf of the number of primary customers that arrive between successive service initiations when X n > 0 (resp.X n = 0).Moreover, we can further assume class dependent service times, i.e., when an orbiting (resp.primary) customer is the one that occupies the server, the pgf of the number of arriving customers during his/her service time equals C o (z) (resp.C p (z)).In such a case C(z) = α1Co(z)+λ1zCp(z) λ1+α1 . Similarly, G(z) = α0Go(z)+λ0zGp(z) λ0+α0 when X n = 0.
Remark 15 Moreover, some very interesting special cases may be deduced from (134).In particular, when α k → ∞, k = 0, 1, then C(z) = C(z), and G(z) = G(z), since α k +λ k z λ k +α k → 1 as α k → ∞.Thus, (134) reduces to the functional equation that corresponds to the standard M/G/1 queue generalization in [1,Section 5].Moreover, one can further assume that one of α k s to tend to infinity, e.g., α 0 → ∞ and a 1 > 0. In such a scenario, the server has the flexibility to treat the orbit queue as a typical queue, when at the beginning of the last service the orbit queue was empty.

An extension to a two-dimensional case: A priority retrial queue
In the following, we go one step further towards a multidimensional case.In particular, we consider the twodimensional discrete-time process {(X 1,n , X 2,n ); n = 0, 1, . ..}, and assume that only the component {X 2,n ; n = 0, 1, . ..} is subject to the autoregressive concept, i.e., we generalize the previous model to incorporate two classes of customers (primary and orbiting customers) and priorities, where orbiting customers are impatient.
Primary customers arrive according to a Poisson process λ 1 and if they find the server busy form a queue waiting to be served.Retrial customers arrive according to a Poisson process λ 2 , and upon finding a busy server join an infinite capacity orbit queue, from where they retry according to the constant retrial policy, i.e., only the first in orbit queue attempts to connect with the server after an exponentially distributed time with rate α.
Let X i,n be the number of customers in queue i (i.e., type i customers) just after the beginning of the nth service, where with i = 1 (i = 2) we refer to the orbit (resp.primary) queue.As usual, the server becomes available to the orbiting customers only when there are no customers at the primary queue upon a service completion.We further assume that orbiting customers become impatient during the service of an orbiting customer, according to the machinery described above.

Conclusion
In this work we investigated the transient and/or the stationary behaviour of various reflected autoregressive processes.These type of processes are described by stochastic recursions where various independence assumptions among the sequences of random variables that are involved there, are lifted and for which, a detailed exact analysis can be also provided.This is accomplished by using Liouville's theorem [11,Theorem 10.52], as well as stating and solving a Wiener-Hopf boundary value problem or an integral equation.Various options for followup research arise.One of them is to concern multivariate extensions of the processes that we introduced.Such vector-valued counterparts are anticipated to be highly challenging.In subsection 8.1 we cope with a simple two-dimensional case, however, the autoregressive parameter was used only in one component.Other possible line of research concerns scaling limits and asymptotics.One also anticipates that, under certain appropriate scalings, a diffusion analysis similar to the one presented in [7] can be applied.