On a class of multiplicative Lindley-type recursions with Markov-modulated dependencies

TL;DR

This paper analyzes Markov-modulated multiplicative Lindley recursions, deriving their stationary distributions and moments, and explores their asymptotic behavior through recursive methods and numerical examples.

Contribution

It introduces new methods to analyze Markov-dependent Lindley-type recursions, deriving explicit Laplace-Stieltjes transforms and recursive formulas for moments.

Findings

Derived the Laplace-Stieltjes transform of the stationary distribution.

Provided recursive formulas for steady-state moments.

Illustrated results with a numerical example.

Abstract

In this paper, we study Markov-modulated dependencies for the multiplicative Lindley's recursion $W_{n+1}=[V_{n}W_{n}+Y_{n}(V_{n})]^{+}$, where $Y_{n}(V_{n})$ may depend on $V_{n}$, and can be written as the difference of two nonnegative random variables that also depend on a common background discrete-time Markov chain $\{Z_{n}\}_{n\in\mathbb{N}}$. Given the state of the background Markov chain, we consider two cases: a) $V_{n}$ equals either 1, or $a\in(0,1)$, or it is negative with certain probabilities, and $Y_{n}(V_{n}):=Y_{n}=S_{n}-A_{n+1}$, where both $A_n$ and $S_n$ have a rational Laplace-Stieltjes transform (LST). b) $V_{n}$ equals $1$ or $-1$ according to certain probabilities, and $Y_{n}(V_{n})$ follow a more general scheme, dependent on $V_{n}$. In both cases, we derive the LST of the stationary transform vector of $\{W_{n}\}_{n\in\mathbb{N}_{0}}$. In the second case, we also provide a recursive approach to obtain the steady-state moments and investigate its asymptotic behavior. A simple numerical example illustrates the theoretical findings.

Equations

$$W_{n+1}=[V_n{W}_n+Y_n(V_n){]}^{+},$$

$$Y_{n}(V_{n})=\left\{\begin{array}[]{ll}S_{n}-A_{n+1},&\text{ when }V_{n}=1,\\ D_{n+1}-C_{n},&\text{ when }V_{n}=-1.\end{array}\right.$$

$$W_{n+1}=[W_n+K_n(V_n)-L_n(V_n){]}^{+},$$

$$\begin{array}[]{rl}K_{n}(V_{n})=&\bar{K}_{n}(V_{n})+\widehat{K}_{n}(V_{n})W_{n},\vskip 5.69054pt\\ L_{n}(V_{n})=&\bar{L}_{n}(V_{n})+\widehat{L}_{n}(V_{n})W_{n},\end{array}$$

$$\begin{array}[]{lr}\bar{K}_{n}(V_{n})=\left\{\begin{array}[]{ll}S_{n},&\text{when }V_{n}=1,\\ D_{n+1},&\text{when }V_{n}=-1,\end{array}\right.,&\bar{L}_{n}(V_{n})=\left\{\begin{array}[]{ll}A_{n+1},&\text{when }V_{n}=1,\\ C_{n},&\text{when }V_{n}=-1.\end{array}\right.\end{array}$$

$$\begin{array}[]{l}P(A_{n+1}\leq x,S_{n}\leq y,Z_{n+1}=j|Z_{n}=i,(C_{r+1},D_{r},Z_{r}),0\leq r\leq n-1)\vskip 5.69054pt\\ =P(A_{n+1}\leq x,S_{n}\leq y,Z_{n+1}=j|Z_{n}=i)\vskip 5.69054pt\\ =p_{i,j}P(A_{n+1}\leq x,S_{n}\leq y|Z_{n}=i,Z_{n+1}=j)=p_{i,j}F_{S,i}(y)G_{A,j}(x),\end{array}$$

$$Y_{n}(V_{n})=\left\{\begin{array}[]{ll}S_{n}-A_{n+1},&\text{ when }V_{n}=1,\\ C_{n+1}-D_{n},&\text{ when }V_{n}=-1.\end{array}\right.$$

$$\begin{array}[]{l}P(D_{n+1}\leq x,C_{n}\leq y,Z_{n+1}=j|Z_{n}=i,(D_{r+1},C_{r},Z_{r}),0\leq r\leq n-1)\vskip 5.69054pt\\ =P(D_{n+1}\leq x,C_{n}\leq y,Z_{n+1}=j|Z_{n}=i)\vskip 5.69054pt\\ =p_{i,j}P(D_{n+1}\leq x,C_{n}\leq y|Z_{n}=i,Z_{n+1}=j)=p_{i,j}F_{C,i}(y)G_{D,j}(x),\end{array}$$

$$W_{n+1}=[V_n{W}_n+S_n-A_{n+1}{]}^{+},$$

$${\Phi }_{B,i}(s):={\int }_0^{\infty }{e}^{-st}d{F}_{S,i}(t):=\frac{N_{B,i}(s)}{D_{B,i}(s)},{\Phi }_{A,j}(s):={\int }_0^{\infty }{e}^{-st}d{G}_{A,i}(t):=\frac{N_{A,j}(s)}{D_{A,j}(s)},$$

$${\Phi }_{W,j}(s)-\sum _{i=1}^N{p}_{i,j}{H}_{i,j}(s)\left[p_1{\Phi }_{W,i}(s)+p_2{\Phi }_{W,i}(as)\right]=\frac{{\sum }_{\omega =0}^{l_j+{\sum }_{k=1}^N{m}_k}{c}_{\omega ,j}{s}^{\omega }}{{\prod }_{k=1}^N{D}_Y^{(k,j)}(s)},$$

$${\tilde{\Phi }}_W(s)=R(s){\tilde{\Phi }}_W(as)+\tilde{V}(s),$$

$$\begin{array}[]{rl}R(s):=&p_{2}(I-p_{1}F(s))^{-1}F(s),\vskip 5.69054pt\\ \tilde{V}(s):=&(I-p_{1}F(s))^{-1}\tilde{v}(s),\end{array}$$

$$\tilde{v}(s):=(\frac{{\sum }_{\omega =0}^{l_1+{\sum }_{k=1}^N{m}_k}{c}_{\omega ,1}{s}^{\omega }}{{\prod }_{k=1}^N{D}_Y^{(k,1)}(s)},\dots ,\frac{{\sum }_{\omega =0}^{l_N+{\sum }_{k=1}^N{m}_k}{c}_{\omega ,N}{s}^{\omega }}{{\prod }_{k=1}^N{D}_Y^{(k,N)}(s)}{)}^T.$$

$$\begin{array}[]{c}\tilde{\Phi}_{W}(s)=\sum_{k=0}^{\infty}\prod_{m=0}^{k-1}R(a^{m}s)\tilde{V}(a^{k}s).\end{array}$$

$$\begin{array}[]{l}\Phi_{W,j}^{n+1}(s)=E\left(e^{-sW_{n+1}}1_{\{Z_{n+1}=j\}}\right)=\sum_{i=1}^{N}P(Z_{n}=i)E\left(e^{-sW_{n+1}}1_{\{Z_{n+1}=j\}}|Z_{n}=i\right)\vskip 5.69054pt\\ =\sum_{i=1}^{N}P(Z_{n}=i)E\left(e^{-s[V_{n}W_{n}+S_{n}-A_{n+1}]^{+}}1_{\{Z_{n+1}=j\}}|Z_{n}=i\right)\vskip 5.69054pt\\ =\sum_{i=1}^{N}P(Z_{n}=i)p_{i,j}E\left(e^{-s[V_{n}W_{n}+S_{n}-A_{n+1}]^{+}}|Z_{n+1}=j,Z_{n}=i\right)\vskip 5.69054pt\\ =\sum_{i=1}^{N}P(Z_{n}=i)p_{i,j}E\left(e^{-s[V_{n}W_{n}+S_{n}-A_{n+1}]}+1-e^{-s[V_{n}W_{n}+S_{n}-A_{n+1}]^{-}}|Z_{n+1}=j,Z_{n}=i\right)\vskip 5.69054pt\\ =\sum_{i=1}^{N}p_{i,j}E\left(e^{-sV_{n}W_{n}}1_{\{Z_{n}=i\}}\right)E\left(e^{-s(S_{n}-A_{n+1})}|Z_{n+1}=j,Z_{n}=i\right)+U^{n}_{j}(s)\vskip 5.69054pt\\ =\sum_{i=1}^{N}p_{i,j}H_{i,j}(s)\left[p_{1}\Phi_{W,i}^{n}(s)+p_{2}\Phi_{W,i}^{n}(as)+p_{3}\int_{-\infty}^{0}\Phi_{W,i}^{n}(sy)P(V_{n}^{-}\in dy)\right]+U^{n}_{j}(s),\end{array}$$

$$U_j^n(s):=\sum _{i=1}^{N}P(Z_n=i)p_{i,j}E\left(1-e^{-s[V_n{W}_n+S_n-A_{n+1}{]}^{-}}|Z_{n+1}=j,Z_n=i\right).$$

$$\begin{array}[]{l}\Phi_{W,j}(s)D_{A,j}(-s)\prod_{k=1}^{N}D_{B,k}(s)\vskip 5.69054pt\\ -\sum_{i=1}^{N}p_{i,j}N_{Y}^{(i,j)}(s)\prod_{\nu\neq i}^{N}D_{B,\nu}(s)\left[p_{1}\Phi_{W,i}(s)+p_{2}\Phi_{W,i}(as)\right]\vskip 5.69054pt\\ =p_{3}\sum_{i=1}^{N}p_{i,j}N_{Y}^{(i,j)}(s)\prod_{\nu\neq i}^{N}D_{B,\nu}(s)\int_{-\infty}^{0}\Phi_{W,i}(sy)P(V^{-}\in dy)\vskip 5.69054pt\\ +D_{A,j}(-s)\prod_{k=1}^{N}D_{B,k}(s)U_{j}(s),\,j=1,\ldots,N,\end{array}$$

$$\begin{array}[]{l}\Phi_{W,j}(s)D_{A,j}(-s)\prod_{k=1}^{N}D_{B,k}(s)-\sum_{i=1}^{N}p_{i,j}N_{Y}^{(i,j)}(s)\prod_{\nu\neq i}^{N}D_{B,\nu}(s)\left[p_{1}\Phi_{W,i}(s)\right.\vskip 5.69054pt\\ \left.+p_{2}\Phi_{W,i}(as)\right]=\sum_{\omega=0}^{l_{j}+\sum_{k=1}^{N}m_{k}}c_{\omega,j}s^{\omega},\end{array}$$

$$\begin{array}[]{l}p_{3}\sum_{i=1}^{N}p_{i,j}N_{Y}^{(i,j)}(s)\prod_{\nu\neq i}^{N}D_{B,\nu}(s)\int_{-\infty}^{0}\Phi_{W,j}(sy)P(V^{-}\in dy)\vskip 5.69054pt\\ +D_{A,j}(-s)\prod_{k=1}^{N}D_{B,k}(s)U_{j}(s)=\sum_{\omega=0}^{l_{j}+\sum_{k=1}^{N}m_{k}}c_{\omega,j}s^{\omega}.\end{array}$$

$$(I-p_{1}F(s)){\tilde{\Phi }}_W(s)=p_{2}F(s){\tilde{\Phi }}_W(as)+\tilde{v}(s),Re(s)\ge 0.$$

$${\tilde{\Phi }}_W(s)=R(s){\tilde{\Phi }}_W(as)+\tilde{V}(s),$$

$${\tilde{\Phi }}_W(s)=\sum _{k=0}^{n-1}\prod _{m=0}^{k-1}R(a^{m}s)\tilde{V}(a^{k}s)+\prod _{m=0}^{n-1}R(a^{m}s){\tilde{\Phi }}_W(a^{n}s).$$

$${\tilde{\Phi }}_W(s)=\sum _{k=0}^{\infty }\prod _{m=0}^{k-1}R(a^{m}s)\tilde{V}(a^{k}s).$$

$$\begin{array}[]{c}\tilde{\zeta}_{r,j}(I-p_{1}F(\delta_{r,j}))=0.\end{array}$$

$$\begin{array}[]{rl}\tilde{\zeta}_{r,j}(I-p_{1}F(\delta_{r,j}))\tilde{\Phi}_{W}(\delta_{r,j})=&\tilde{\zeta}_{r,j}[p_{2}F(\delta_{r,j})\tilde{\Phi}_{W}(a\delta_{r,j})+\tilde{v}(\delta_{r,j})]=0,\end{array}$$

$$p_2{\tilde{\zeta }}_{r,j}F({\delta }_{r,j}){\tilde{\Phi }}_W(a{\delta }_{r,j})=-\tilde{\zeta }\tilde{v}({\delta }_{r,j}).$$

$$D_{A,j}(-s)\prod _{k=1}^N{D}_{B,k}(s)=\prod _{r=1}^{l_j}(-s-s_{r,j})\prod _{k=1}^N\prod _{f=1}^{m_k}(t-t_{f,k}).$$

$$\begin{array}[]{l}p_{3}\sum_{i=1}^{N}p_{i,j}N_{Y}^{(i,j)}(t_{f,k})\prod_{\nu\neq i}^{N}D_{B,\nu}(t_{f,k})\int_{-\infty}^{0}\left[\sum_{v=0}^{\infty}\prod_{m=0}^{v-1}R(a^{m}t_{f,k}y)\tilde{V}(a^{v}t_{f,k}y)\right]\vskip 5.69054pt\\ \times P(V^{-}\in dy)=\sum_{\omega=0}^{l_{j}+\sum_{k=1}^{N}m_{k}}c_{\omega,j}(t_{f,k})^{\omega},\,j=1,\ldots,N.\end{array}$$

Full text

Keywords: Workload; Laplace-Stieltjes transform; Recursion; Markov-modulation

1 Introduction

In this work, we focus on the stationary analysis of the Markov-modulated multiplicative version of the Lindley-type recursion

[TABLE]

with $[x]^{+}:=max\{x,0\}$. We assume that $Y_{n}(V_{n})$ can be written as a difference of two nonnegative random variables that may depend on $V_{n}$. Stochastic recursions have received much attention in the applied probability literature, and especially those of autoregressive type, due to their wide applicability across various scientific domains including biology, finance, and engineering, e.g., [12, 13, 20, 16, 24].

The vast majority of the related works focus on the investigation of the stability conditions, limit theorems and questions related to queuing applications [30], whereas our primary focus lies on the derivation of results for the stationary distribution of the process under investigation. In [5, 4, 21], the authors studied similar linear stochastic recursions (with applications in branching processes, and queueing systems), but their focus was only on deriving the first two (stationary, and transient) moments, and not on deriving the distribution of their main argument. Furthermore, most of the works on this field refers to the scalar case. Our focus is on the derivation of the stationary distribution (in terms of LST) of vector-valued stochastic recursions of the type in (1), which in turn result in vector-valued functional equations.

In this work, we focus on two models. In Model I, the $V_{n}$ are either equal to one, or equal to a positive constant $a\in(0,1)$, or negative, according to certain probabilities. We further assume that $Y_{n}(V_{n})=S_{n}-A_{n+1}$ (i.e., independent of $V_{n}$), and we demand that both the positive and the negative parts of the $Y_{n}$ have a rational Laplace–Stieltjes transform (LST). Moreover, $S_{n}$, $A_{n+1}$ are regulated by a finite state space irreducible Markov chain, thus, given the state of the background Markov chain, $S_{n}$, $A_{n+1}$ are conditionally independent.In Model II, the $V_{n}$ are equal either to one, or to minus one, and $Y_{n}(V_{n})$ have a more general form that is dependent on $V_{n}$.

1.1 Related work

In [10], the authors considered the (scalar) recursion $W_{n+1}=[aW_{n}+S_{n}-A_{n}]^{+}$, where $\{S_{n}-A_{n}\}_{n\in\mathbb{N}_{0}}$ forms a sequence of independent and identically distributed (i.i.d.) random variables and $a\in(0,1)$. The authors investigated both the transient and stationary behaviour. Note that in such a case, $W_{n}$ could be interpreted as the workload in a queueing system just before the $n$th arrival, which adds $S_{n}$ work, and makes obsolete a fixed fraction $1-a$ of the already present work. 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$ was investigated in [28].

In [8], the authors provided transient and the stationary characteristics for the scalar autoregressive process described by the recursion $W_{n+1}=[V_{n}W_{n}+S_{n}-A_{n}]^{+}$, where $V_{n}$ can be positive or negative, while in [7], the authors considered the case where $V_{n}W_{n}$ was replaced by $S(W_{n})$, where $\{S(t)\}$ is a Levy subordinator. Motivated by applications that arise in queueing and insurance risk models, the authors in [9] considered Lindley-type recursions where $\{B_{n}\}_{n\in\mathbb{N}_{0}}$, $\{A_{n}\}_{n\in\mathbb{N}_{0}}$ obey a semi-linear dependence. The seminal work in [1] provided a method to study functional equations that arise in a wide range of queueing, autoregressive and branching processes, and where (scalar) Lindley-type recursions arise. A generalized version of the model in [8], was investigated in [22]. Quite recently, in [19] the authors generalized the work in [10] by considering, among others, non-trivial dependence structures among $\{S_{n}\}_{n\in\mathbb{N}}$, $\{A_{n}\}_{n\in\mathbb{N}}$. The model in [11] is related to our work (in Model II) but corresponds to the scalar case.

Quite recently, in [18], the author investigated vector-valued reflected autoregressive processes that are described by vector-valued Lindley type recursions of the form $W_{n+1}=[aW_{n}+S_{n}-A_{n+1}]^{+}$, where the sequences $\{S_{n}\}_{n\in\mathbb{N}_{0}}$, $\{A_{n}\}_{n\in\mathbb{N}_{0}}$ are governed by an irreducible finite state background Markovian process, thus, considered the Markov-dependent version of the process analysed in [10]. Note that the specific case of $a=1$ in [18] corresponds to the waiting time in a single server queue with Markov-dependent interarrival and service times studied in [2]; see also [3] for a Markov-dependent structure related in insurance mathematics. The case where $a=-1$ was investigated in [29] in the context of carousel models. In [18] the focus was on the case where $a\in(0,1)$. Moreover, the author also discussed the case where $\{(S_{n},A_{n+1})\}_{n\in\mathbb{N}_{0}}$ form a sequence of i.i.d. random vectors with a distribution function defined by the Farlie-Gumbel-Morgenstern copula, as well as the case where $\{(S_{n},A_{n+1})\}_{n\in\mathbb{N}_{0}}$ have a bivariate matrix-exponential distribution. In both cases, the joint distribution is also dependent on the state of the background discrete time Markov chain.

1.2 Contribution

A major contribution of the present paper, regarding Model I, is that we derive and solve a system of Wiener–Hopf equations, which allows us to study the stationary behavior of the $\{W_{n}\}_{n\in\mathbb{N}_{0}}$ process. To the best of the author’s knowledge, in Model I we studied for the first time a recursion of the form $W_{n+1}=[V_{n}W_{n}+Y_{n}(V_{n})]^{+}$, where $Y_{n}(V_{n})=Y_{n}=S_{n}-A_{n+1}$, $P(V_{n}=1)=p_{1}$, $P(V_{n}=a)=p_{2}$, where $a\in(0,1)$, $P(V_{n}<0)=p_{3}=1-p_{1}-p_{2}$, under a semi-Markov dependent framework introduced in [2] (see Section 3). We focus on the most interesting case where $0<p_{i}<1$, $i=1,2,3$, thus, we provide a significant extension to the existing results by allowing the $V_{n}$ to take arbitrary distributions on $(-\infty,0]$. Here, we demand that both the positive and the negative parts of the $Y_{n}$ have a rational LST. Moreover, contrary to the work in [8] where the authors studied a scalar version (non-modulated), we introduce Markov-modulated dependencies, which lead to a vector-valued fixed-point functional equation that will be solved recursively, with some additional technical requirements. By studying (1), we significantly extend the analysis of the Lindley recursion as well as the analysis of the Markov-modulated stochastic recursion in [18], where vector-valued recursions of the above form were studied for the case $P(V_{n}=a)=1$, $a\in(0,1)$. By letting the random variable $V_{n}$ to be negative we introduce to the model specific mathematical intricacies, but still, we are able to get explicit results in terms of Laplace transforms. In contrast to the scalar case analyzed in [8], our vector-valued framework results in a system of Wiener-Hopf equations (instead of a single equation in [8]) that we succeed to solve and prove the convergence of its solution. We also discuss the stability condition. Moreover, contrary to [21], where the authors focused on the derivation of the first two moments, in our work our primary focus is on the stationary distribution of the process under study.

Moreover, we also analyzed a special case that results in a different solution approach. In particular, we consider the case where $p_{3}=1$ (i.e., $p_{1}=p_{2}=0$), and we drop the assumption that $A_{n}$, given the state of the background Markov chain, has a rational LST. With this special case, we generalize the corresponding model in [8, Section 4] from the scalar case to the vector case, thus, filling an important gap in the literature, since we are now dealing with a vector-valued functional equation. 2. 2.

Our second contribution refers to the analysis of a Markov-modulated Lindley-type recursion in Model II, for which the monotonicity assumption does not hold. The scalar Lindley recursion is well studied in the applied probability literature and refers to the class of discrete time Markov chains described by the recursion $W_{n+1}=F(W_{n},X_{n})$. An important assumption that is often made in the related literature is that the function $F(w,x)$ is non-decreasing in its main argument $w$ [6]. In particular, the recursion we focus on does not satisfy that usual assumption, which in turn, does not allow us to exploit the known duality results between Markov processes and monotone increasing continuous processes; [6, Section IX.4]. More precisely, we go one step further and study in detail the Markov-modulated version of that recursion that does not satisfy the monotonicity property, thus, filling an important gap in the related literature. Moreover, although we incorporate Markov-modulation to capture fluctuations in interarrival and service times, we are still able to solve explicitly the vector-valued functional equation in terms of Laplace transforms. Contrary to the scalar case in [11], where the authors by solving a Wiener-Hopf boundary value problem derived the LST of the waiting time distribution, we are able to solve a vector-valued functional equation by using simpler arguments, and without using the Wiener-Hopf method. This is quite promising since our Markov-modulated framework we use, allows for solving even more complicated problems without increasing significantly the mathematical difficulties.

More precisely, in Section 4 we study the Markov-dependent multiplicative Lindley-type recursion

[TABLE]

where for every $n$ the random variable $V_{n}$ is equal to plus or minus one according to $P(V_{n}=1)=p$, $P(V_{n}=-1)=q:=1-p$, $p\in[0,1]$. Moreover, for each $n$ the random variable $Y_{n}(V_{n})$ is written as the difference of two random variables as shown below:

[TABLE]

On top of that, we assume that the sequences $\{A_{n}\}_{n\in\mathbb{N}_{0}}$, $\{B_{n}\}_{n\in\mathbb{N}}$ (resp. $\{C_{n}\}_{n\in\mathbb{N}_{0}}$, $\{D_{n}\}_{n\in\mathbb{N}}$) are assumed to be conditional independent sequences of i.i.d. non-negative random variables, given the state of an irreducible background Markov chain with finite state space (i.e., Markov-modulated dependencies). Our main goal is to derive the Laplace-Stieltjes transform (LST) of the steady-state distribution of $\{W_{n},n=1,2,...\}$, when it exists.

To the author’s best knowledge, it is the first time in the related literature that such a Markov-modulated multiplicative version of Lindley recursion is studied. Contrary to the case in [11], where the authors studied a scalar version of (2), in our work we consider the Markov-modulated version, by further allowing $Y_{n}$ to depend on $V_{n}$. For such a model, we derive the transform vector of the LSTs of the steady-state distribution of the waiting time. We further provided a recursive approach to obtain the steady-state moments of the waiting time distribution, as well as to investigate the asymptotic behaviour of the steady-state waiting distribution for light-tailed service time distributions.

Despite our methodological contribution, our work has also strong connections with respect to the performance analysis of service systems. In particular, this state-dependent queuing processes, as discussed below (see also [11] for the scalar case), are connected to LaPalice queues [23], in which customers are scheduled to be served such that the period between two consecutively scheduled customers is greater than or equal to the service time of the first customer. The Markov-modulated framework we incorporate, provide an even more realistic aspect of the model, since in many practical situations, the customer arrivals and the service mechanism of a queueing system are both influenced by some factors external or internal to the system [25]. In financial terms, $W_{n}$ may represent an inventory in time period $n$ (e.g. cash), $V_{n}$ may represent a multiplicative, possibly random, decay or growth factor between times $n$ and $n+1$ (e.g. interest rate) and $Y_{n}(V_{n})$ may represent a quantity that is added or subtracted between times $n$ and $n+1$ (e.g. deposit minus withdrawal).

When $V_{n}=1$ (i.e., $p=1$) we have the MAP/G/1 queue [2], and when $V_{n}=-1$ (i.e., $p=0$, equivalently $q=1$) we have a modulated alternating service model with two service points that arise in the performance analysis of warehouse systems; see [26, 29]. Studying a recursion that contains the Markov-modulated version of both Lindley’s classical recursion [2] and the recursion in [26, 28] as special cases it is of interest in its own right. For $p\in(0,1)$, our model describes a Markov-modulated FCFS queue, in which the service times and the interarrival times depend linearly and randomly on the waiting times, as well as on the value of $V_{n}$.

We now discuss this special queueing model in more detail. Consider an extension of a modulated G/G/1 queue in which the service times and the interarrival times depend linearly and randomly on the waiting times. In particular, the model is specified by a stationary and ergodic sequence of non-negative random variables $\{(\bar{K}_{n}(V_{n}),\bar{L}_{n}(V_{n}),\widehat{K}_{n}(V_{n}),\widehat{L}_{n}(V_{n}),V_{n}\}$, $n\geq 0$. Given the state of an irreducible finite state Markov chain, the sequence $\{W_{n};n\geq 0\}$ is defined recursively by

[TABLE]

where

[TABLE]

with

[TABLE]

Note that $\bar{K}_{n}(V_{n})$ can be considered as the nominal service time of the $n$ customer, while $\bar{L}_{n}(V_{n})$ can be considered as the interarrival time between customers $n$ and $n+1$ customer, and are both dependent on $V_{n}$, as well as on the state of the background Markov chain.

Then, we have a modulated multiplicative Lindley recursion with $V_{n}=1+\widehat{K}_{n}(V_{n})-\widehat{L}_{n}(V_{n})$. When $V_{n}=1$ (with probability $p$), then $\widehat{K}_{n}(1)=\widehat{L}_{n}(1)$, while when $V_{n}=-1$ (with probability $q=1-p$), $\widehat{L}_{n}(-1)=2+\widehat{K}_{n}(-1)$. This state-dependent queuing processes, as also mentioned in [11], are connected to a class of queueing systems [23], in which customers are scheduled such that the period between two consecutively scheduled customers is greater than or equal to the service time of the first customer.

1.3 Structure of the paper

The paper is organized as follows. In Section 2 we provide the description of the two models and some preliminaries. Sections 3, 4 are devoted to the stationary analysis of Models I, II, respectively. In Section 5, we provide a simple numerical example to illustrate the theoretical findings. Some concluding remarks and suggestions for future research are presented in Section 6.

2 Model description and preliminaries

In this work we focus on the stationary analysis of a Markov-dependent version of the recursion (1), where $\{Y_{n}(V_{n})\}_{n\in\mathbb{N}_{0}}$ is a sequence of i.i.d random variables that is written as a difference of two random variables. Their distributions are regulated by an irreducible discrete-time Markov chain $\{Z_{n}\}_{n\in\mathbb{N}_{0}}$ with state space $E=\{1,2,...,N\}$, one-step transition probability matrix $P:=(p_{i,j})_{i,j\in E}$, and stationary distribution $\tilde{\pi}:=(\pi_{1},\ldots,\pi_{N})^{T}$.

Let $\Phi_{W,i}^{n}(s):=E(e^{-sW_{n}}1_{\{Z_{n}=i\}})$, $i=1,\ldots,N,$ $Re(s)\geq 0$, $n\geq 0$, and assuming the limit exists, define $\Phi_{W,i}(s)=\lim_{n\to\infty}\Phi_{W,i}^{n}(s)$, $i=1,\ldots,N$. Let also $\tilde{\Phi}_{W}(s):=(\Phi_{W,1}(s),\ldots,\Phi_{W,N}(s))^{T}$.

In this paper, we discuss the following two variants of the model in (1):

Model I: $P(V_{n}=1)=p_{1}$, $P(V_{n}=a)=p_{2}$, $a\in(0,1)$, $P(V_{n}<0)=p_{3}:=1-p_{1}-p_{2}$, $Y_{n}(V_{n}):=Y_{n}=S_{n}-A_{n+1}$. We further assume that for $n\geq 0$, $x,y\geq 0$, $i,j=1,\ldots,N$:

[TABLE]

where $F_{S,i}(.)$, $G_{A,j}(.)$ denote the distribution functions of $S_{n}$, $A_{n+1}$, given $Z_{n}=i$, $Z_{n+1}=j$, respectively. Note that $A_{n+1}$, $S_{n}$, $Z_{n+1}$ are independent of the past given $Z_{n}$, and $A_{n+1}$, $S_{n}$ are conditionally independent given $Z_{n}$, $Z_{n+1}$. We assume that the sequences $\{A_{n}\}_{n\in\mathbb{N}}$ and $\{S_{n}\}_{n\in\mathbb{N}_{0}}$ are autocorrelated as well as cross-correlated. 2. 2.

Model II: $P(V_{n}=1)=p$, $P(V_{n}=-1)=q:=1-p$, and

[TABLE]

Note that when $V_{n}=1$, the dynamics among $\{A_{n}\}_{n\in\mathbb{N}}$, $\{S_{n}\}_{n\in\mathbb{N}_{0}}$ are governed by (4). When $V_{n}=-1$, then,

[TABLE]

where $F_{C,i}(.)$, $G_{D,j}(.)$ denote the distribution functions of $C_{n}$, $D_{n+1}$, given $Z_{n}=i$, $Z_{n+1}=j$, respectively. Note that $D_{n+1}$, $C_{n}$, $Z_{n+1}$ are independent of the past given $Z_{n}$, and $D_{n+1}$, $C_{n}$ are conditionally independent given $Z_{n}$, $Z_{n+1}$. We assume that the sequences $\{D_{n}\}_{n\in\mathbb{N}}$ and $\{C_{n}\}_{n\in\mathbb{N}_{0}}$ are autocorrelated as well as cross-correlated.

3 Model I: The Markov-dependent mixed case

Starting from the recursion

[TABLE]

and under the Markov-dependent dynamics in (4), we assume that $V_{n}=1$ with probability $p_{1}$, $V_{n}=a\in(0,1)$ with probability $p_{2}$, and $V_{n}<0$ with probability $p_{3}=1-p_{1}-p_{2}$. We consider the case where $0<p_{i}<1$, $i=1,2,3,$ $p_{1}+p_{2}<1$. If $p_{1}=1$, the model is reduced to the one in [2], while if $p_{2}=1$ is reduced to the one in [18], and if $p_{3}=1$ is reduced to the one in [29].

3.1 On the stability condition

We denote with $X$ a generic random variable distributed as $X_{0}$, and always assume that $E(|V|)<\infty$, $E(|Y|)<\infty$. For the scalar case, in [30, Theorem 1] it is shown that if $P(V<0)>0$ and $P(Y\leq 0)>0$, then, $W_{n}$ tends to a proper limit as $n\to\infty$. Note that for our case, $V$ does not depend on the state of the background discrete-time Markov chain, although $P(V<0)=p_{3}>0$. The other requirement is always satisfied since, in the simple case where $A|Z=j\sim exp(\lambda_{j})$, and $S|Z=j$ follow a general distribution with distribution function $B_{j}(.)$, and LST $\beta_{j}^{*}(.)$, $P(Y\leq 0)=\sum_{i=1}^{N}\sum_{j=1}^{N}\pi_{i}p_{i,j}\beta_{i}^{*}(\lambda_{j})>0$, thus, from the regenerative nature of $W_{n}$, the limit $W$ is unique.

3.2 Stationary analysis

We further assume that $S|Z=i$ and $A|Z=j$ have a rational LST, i.e.,

[TABLE]

where $\Phi_{B,i}\in\mathbb{Q}[t_{1,i},\ldots,t_{m_{i},i}]$, $i=1,\ldots,N$, $\Phi_{A,j}\in\mathbb{Q}[s_{1,j},\ldots,s_{l_{j},j}]$, $j=1,\ldots,N,$ with $Re(s_{r,j})<0$, $r=1,\ldots,l_{j}$, $j=1,\ldots,N$ and $Re(t_{k,i})<0$, $k=1,\ldots,m_{i}$, $i=1,\ldots,N$. The rationality assumptions are natural in the related literature for the non-modulated G/G/1 queue. While in principle the waiting-time distribution in the general standard G/G/1 queue can be obtained via a Wiener–Hopf decomposition [15, Chapter II.5], the solution is a rather implicit one, unless one makes rationality assumptions on either the interarrival or the service-time LST [17]. Moreover, it is well-known that the distribution of any nonnegative random variable can be approximated arbitrarily closely by the distribution of a random variable with a rational LST [6, Ch. III], so that a restriction to random variables with rational LST leads to just a minor loss of generality.

In such a setting, denote by $H_{i,j}(s):=E(e^{-sY_{n}}|Z_{n}=i,Z_{n+1}=j)=E(e^{-s(S_{n}-A_{n+1})}|Z_{n}=i,Z_{n+1}=j)=\frac{N_{Y}^{(i,j)}(s)}{D_{Y}^{(i,j)}(s)}=\frac{N_{B,i}(s)N_{A,j}(-s)}{D_{B,i}(s)D_{A,j}(-s)}.$ Denote further the $N\times N$ matrix $H(s):=(H_{i,j}(s))_{i,j=1,\ldots,N}$. Then, we have the following result.

Theorem 1

For $s\in\Omega:=\{s:Re(s)\geq 0,det(I-p_{1}F(s))\neq 0\}$, the transforms $\Phi_{W,j}(s)$, $j=1,\ldots,N$ satisfy the system

[TABLE]

where $c_{0,j}=\pi_{j}p_{3}\prod_{k=1}^{N}D_{Y}^{(k,j)}(0)$, and $c_{\omega,j}$, $\omega=1,\ldots,l_{j}+\sum_{k=1}^{N}m_{k}$, $j=1,\ldots,N$ can be determined from the linear systems (18), (19) given below.

Equivalently, in matrix notation, the transform vector $\tilde{\Phi}_{W}(s)$ satisfies

[TABLE]

where

[TABLE]

and $F(s):=P^{T}\circ H(s)$, with ’$\circ$’ denotes the Hadamard product and

[TABLE]

Then, the solution to (8) is

[TABLE]

Proof. Using the recursion (6), and the identity $e^{-s[x]^{+}}+e^{-s[x]^{-}}=e^{-sx}+1$ (where $[x]^{+}:=max\{x,0\}$, $[x]^{-}:=min\{x,0\}$), we have for $Re(s)=0$, the following equation for the transforms $\Phi_{W,j}^{n}(s)$, $j=1,\ldots,N$:

[TABLE]

where $V_{n}^{-}=\left(V_{n}|V_{n}<0\right)$, and

[TABLE]

Letting $n\to\infty$, substituting $H_{i,j}(s)$ and rearranging we come up with the following functional equation for $Re(s)=0$:

[TABLE]

where $U_{j}(s)=\lim_{n\to\infty}U^{n}_{j}(s)$. Now the following are true:

the left-hand side of (10) is analytic in $Re(s)>0$ and continuous in $Re(s)\geq 0$, 2. 2.

the right-hand side of (10) is analytic in $Re(s)<0$ and continuous in $Re(s)\leq 0$, 3. 3.

for large $s$, and $j=1,\ldots,N,$ both sides are $\mathcal{O}(s^{l_{j}+\sum_{k=1}^{N}m_{k}})$ in their respective half-planes.

Both sides in (10) are well defined at the boundary $Re(s)=0$, and the determination of the unknown functions $\Phi_{W,j}(s)$ and $U_{j}(s)$ under the conditions 1., 2., and 3. is reduced to the solution of a Wiener–Hopf boundary value problem; see [14]. By introducing a function $K(s)$ that is equal to the left-hand side of (10) for $Re(s)\geq 0$ and to the right-hand side of (10) for $Re(s)\leq 0$, we have a function that is analytic in the whole $s$-plane, and that for large $s$ is $\mathcal{O}(s^{l_{j}+\sum_{k=1}^{N}m_{k}})$. Then, Liouville’s theorem [27, p. 85] states that in their respective half-planes both sides of (10) are equal to the same $(s^{l_{j}+\sum_{k=1}^{N}m_{k}})-$th degree polynomial in $s$. Thus, for $Re(s)\geq 0$,

[TABLE]

and, for $Re(s)\leq 0$,

[TABLE]

For $s=0$, in either (11), or (12) yields $c_{0,j}=\pi_{j}p_{3}D_{A,j}(0)\prod_{k=1}^{N}D_{B,k}(0)$. Note that to fully specify $\tilde{\Phi}_{W}(s)$, we need to derive $c_{\omega,j}$, $\omega=1,2,\ldots,l_{j}+\sum_{k=1}^{N}m_{k}$, $j=1,\ldots,N$, i.e., we need to derive a total $\sum_{j=1}^{N}l_{j}+N\sum_{k=1}^{N}m_{k}$ equations for the unknown terms $c_{\omega,j}$.

This task can be accomplished in the following steps:

First note that (11) can be rewritten as (7). In matrix notation, (11) (or equivalently, (7)) results in

[TABLE]

For $s\in\Omega=\{s:Re(s)\geq 0,det(I-p_{1}F(s))\neq 0\}$,

[TABLE]

where $R(s)=(I-p_{1}F(s))^{-1}p_{2}F(s)$, $\tilde{V}(s)=(I-p_{1}F(s))^{-1}\tilde{v}(s)$. Iterating (14) $n$ times yield

[TABLE]

As $n\to\infty$, $\tilde{\Phi}_{W}(a^{n}s)\to\tilde{\pi}$ and it is readily seen that as $m\to\infty$, $R(a^{m}s)\to p_{2}(I-p_{1}P^{T})^{-1}P^{T}$. For $p_{1}+p_{2}<1$, note also that $||p_{2}(I-p_{1}P^{T})^{-1}P^{T}||\leq r<1$ (where $||C||=max_{i,j}|C_{i,j}|$ is a matrix norm of $C\in\mathbb{C}^{N\times N}$). Therefore, $p_{2}(I-p_{1}P^{T})^{-1}P^{T}$ is a convergent matrix, i.e., as $m\to\infty$, $[p_{2}(I-p_{1}P^{T})^{-1}P^{T}]^{m}\to O$ (where $O$ denotes the zero matrix). Therefore, letting $n\to\infty,$

[TABLE] 2. 2.

Following an idea in [17, Appendix 2] (see also [2], and Appendix A), we can show that $det(I-p_{1}F(s))=0$ has exactly $\sum_{j=1}^{N}l_{j}$ roots, say $\delta_{r,j}$, $r=1,\ldots,l_{j}$, $j=1,\ldots,N$, such that $Re(\delta_{r,j})>0$. Therefore, there exist nonzero row vectors $\tilde{\zeta}_{r,j}$, $r=1,\ldots,l_{j}$, $j=1,\ldots,N$, such that

[TABLE]

Substituting $s=\delta_{r,j}$, $r=1,\ldots,l_{j}$, $j=1,\ldots,N$ in (13), pre-multiply the resulting equation with $\tilde{\zeta}_{r,j}$, and having in mind (17) results in

[TABLE]

or equivalently,

[TABLE]

Substituting the right-hand side of (16) for $s=a\delta_{r,j}$ in (18) gives a system of $\sum_{j=1}^{N}l_{j}$ linear equations. 3. 3.

The remaining $N\sum_{k=1}^{N}m_{k}$ are obtained as follows: Note first that

[TABLE]

Setting $t=t_{f,k}$, $f=1,\ldots,m_{k}$, $k=1,\ldots,N,$ in (12), and substituting (16) in the resulting equation yields

[TABLE]

By solving (18), (19), we obtain $c_{\omega,j}$, $\omega=1,2,\ldots,l_{j}+\sum_{k=1}^{N}m_{k}$, $j=1,\ldots,N$, so that $\tilde{\Phi}_{W}(s)$ is fully specified by (16). This finishes the proof.

Once $\tilde{\Phi}_{W}(s)$ is fully specified, we can obtain the steady-state moments. Let $m_{i}:=\lim_{n\to\infty}E(W1_{\{Z_{n}=i\}})$, $i=1,\ldots,N$, and $\tilde{m}=(m_{1},\ldots,m_{N})^{T}$. Differentiating (13) with respect to $s$ and letting $s\to 0$ we obtain

[TABLE]

where $\Phi=P^{T}\circ H^{\prime}(0)$, where $"\prime"$ denotes the first derivative.

3.3 A special case

We briefly describe a special case of Model I by assuming $p_{3}=1$, i.e., $V<0$ a.s.. Moreover, we drop the assumption of rationality of $\Phi_{A}(.)$. We prefer to give a separate analysis of this special case, to make the reader familiar with the specific mathematical intricacies due to $V_{n}$ being negative and $V_{n}$ being a positive constant (see the model in [18]). Moreover, the additional assumption that $S_{n,i}$, $A_{n,i}$ have a rational LST plays an important role in the analysis of Model I, excluding the possibility of simply obtaining the results for this special case directly from those for Model I. Then, following similar arguments as in Theorem 1 we obtain for $Re(s)=0$,

[TABLE]

Let’s first discuss the stability properties of the system. We assume that there are $i$ and $j$ such that $P(S_{n}<A_{n+1},Z_{n+1}=l|Z_{n}=i)>0$. This implies that the Markov chain $\{(W_{n},Z_{n}),n\in\mathbb{N}_{0}\}$ is stable; see also [29] for a similar argument. So, since $\{Z_{n}\}$ is aperiodic, the limiting distribution of $\{(W_{n},Z_{n}),n\in\mathbb{N}_{0}\}$ exists, and $Re(s)>0$, $n>1$, and $j=1,2,\ldots,N$, $\Phi_{W,j}^{n}(s)\to\Phi_{W,j}(s)$, as $n\to\infty$. Therefore, for $n\to\infty$, $Re(s)=0$,

[TABLE]

Then,

the left-hand side of (21) is analytic in $Re(s)>0$ and continuous in $Re(s)\geq 0$, 2. 2.

the right-hand side of (21) is analytic in $Re(s)<0$ and continuous in $Re(s)\leq 0$, 3. 3.

for large $s$, and $j=1,\ldots,N,$ both sides are $\mathcal{O}(s^{\sum_{k=1}^{N}m_{k}})$ in their respective half-planes.

Therefore, for $Re(s)\geq 0$,

[TABLE]

and for $Re(s)\leq 0$,

[TABLE]

For $s=0$, (22), or (23) yields $c_{0,j}=\pi_{j}\prod_{k=1}^{N}D_{B,k}(0).$ We still need $N\sum_{k=1}^{N}m_{k}$ unknown terms. These terms are obtained as follows:

Substituting $t=t_{f,k}$, $f=1,\ldots,m_{k}$, $k=1,\ldots,N,$ in (23), and substituting (22) into the integral (recall that $Re(t_{f,k}<0)$, $\forall f,k$) the resulting equation yields

[TABLE]

The linear system of equations (24) will provide the unknowns $c_{\omega,j}$, $\omega=1,\ldots,\sum_{k=1}^{N}m_{k}$, $j=1,\ldots,N$, so that for $Re(s)\geq 0$,

[TABLE]

4 Model II: On a Markov-modulated queue with service and interarrival times depending on waiting times

We now focus on Markov-modulated recursions of the form

[TABLE]

where $P(V_{n}=1)=p$, $P(V_{n}=-1)=q:=1-p$, and

[TABLE]

In this work, we assume that $A_{n}|Z_{n}=j\sim exp(\lambda_{j})$, $D_{n}|Z_{n}=j\sim exp(\mu_{j})$, while $S_{n}|Z_{n}=j$, $C_{n}|Z_{n}=j$, are arbitrarily distributed with distribution functions $B_{j}(.)$, and $C_{j}(.)$, $j=1,\ldots,N,$ respectively. Moreover, denote their LSTs with $\beta_{j}^{*}(s)=\int_{0}^{\infty}e^{-st}dB_{j}(t)$, $c_{j}^{*}(s)=\int_{0}^{\infty}e^{-st}dC_{j}(t)$. Assume also that $\gamma_{j}=\int_{0}^{\infty}xdB_{j}(x)$, $\delta_{j}=\int_{0}^{\infty}xdC_{j}(x)$, $j=1,\ldots,N$.

Our focus here is on the case where $p\in(0,1)$, so that the service times and the interarrival times depend linearly and randomly on the waiting times, as well as on the value of $V_{n}$.

4.1 On the stability condition

Following ideas from [30, 11], note that from (26), by replacing $V_{n}$ with $max\{0,V_{n}\}$, and $Y_{n}(V_{n})$ with $max\{0,Y_{n}(V_{n})\}$, then, the resulting waiting times are least as large as those given in (26). Then, exploiting ideas from [30], if $W_{n}$ satisfies (26), then with probability 1, $W_{n}\leq Q_{n}$ for all $n$, where

[TABLE]

and $Q_{0}=W_{0}\geq 0$. Thus, if $W_{n}$ satisfies (26), $Q_{n}$ satisfies (27), and $Q_{n}$ converges to the proper limit $Q$, then $\{W_{n}\}_{n\in\mathbb{N}}$ is tight and $P(W>x)\leq P(Q>x)$ for all $x$, where $W$ is the limit in distribution of any convergent subsequence of $\{W_{n}\}_{n\in\mathbb{N}}$. With that in mind, as well as [12, Theorem 1], implies that $Q_{n}$ satisfying (27) converges to a proper limit if $P(max\{0,V_{n}\}=0)=P(V_{n}\leq 0)=q>0$. Then, looking at [30, Theorem 1], the sequence $\{W_{n}\}_{n\in\mathbb{N}}$ is tight for all $\rho=\frac{E(\bar{K}_{0}(V_{0}))}{E(\bar{L}_{0}(V_{0}))}=\frac{\sum_{i=1}^{N}\pi_{i}(p\gamma_{i}+q\mu_{i}^{-1})}{\sum_{i=1}^{N}\pi_{i}(p\lambda_{i}^{-1}+q\delta_{i})}$ and $W_{0}$. If, in addition, $0\leq p<1$ and $\{(V_{n},Y_{n}(V_{n}))\}_{n\in\mathbb{N}_{0}}$ given the state of the background discrete time Markov chain is a sequence of independent vectors with $P(V_{0}\leq 0,Y_{0}(V_{0}))>0$, then the events $\{W_{n}=0\}$ are regeneration points with finite mean time and $\{W_{n}\}_{n\in\mathbb{N}_{0}}$ converges in distribution to a proper limit $W$ as $n\to\infty$ for all $\rho$ and $W_{0}$. Note that when $p=1$, we need to have $\rho<1$; [2].

4.2 Stationary analysis

Recall that $\Phi_{W,i}^{n}(s)=E(e^{-sW_{n}}1_{\{Z_{n}=i\}})$, $Re(s)\geq 0$, $n\geq 0$, $i=1,\ldots,N$, and assuming the limit exists, let $\Phi_{W,i}(s)=\lim_{n\to\infty}\Phi_{W,i}^{n}(s)$. Some additional notation: $\Lambda=diag(\lambda_{1},\ldots,\lambda_{N})$, $M=diag(\mu_{1},\ldots,\mu_{N})$, and

[TABLE]

Theorem 2

The transform row vector $\tilde{\Phi}_{W}^{T}(s)$ satisfies

[TABLE]

where, $\tilde{\Phi}_{W}^{T}(0)=\tilde{\pi}^{T}$,

[TABLE]

Proof. The transform $\Phi_{W,j}^{n+1}(s):=E(e^{-sW_{n+1}}1_{\{Z_{n+1}=j\}})$ satisfies

[TABLE]

Using standard arguments and following [2, 29] we have

[TABLE]

Substituting back in (29), and letting $n\to\infty$ so that $\lim_{n\to\infty}\Phi_{W,i}^{n}(s)=\Phi_{W,i}(s)$, we obtain after some algebra

[TABLE]

In matrix form equation (30) yields (28). Moreover, $\Phi_{W,j}^{n+1}(0)=E(1_{\{Z_{n+1}=j\}})$, and as $n\to\infty$, $\Phi_{W,j}^{n}(0)\to\pi_{j}$, thus, $\tilde{\Phi}^{T}_{W}(0)=\tilde{\pi}^{T}$. This completes the proof.

Note that in order to fully specify $\tilde{\Phi}^{T}_{W}(s)$ we first need to obtain $2N^{2}$ unknown terms, namely $\Phi_{W,i}(\lambda_{j})$, $\Phi_{W,i}(\mu_{j})$, $i,j=1,\ldots,N$, which are involved in the vectors $\tilde{k}_{1}$, $\tilde{k}_{2}$. This task is accomplished by performing the following steps.

First, the following theorem gives the number and placement of the solutions of $det(H(s)+sI-\Lambda)=0$.

Theorem 3

The equation

[TABLE]

has exactly $N$ zeros, say $s_{1},\ldots,s_{N}$, with $Re(s_{i})>0$, $i=1,\ldots,N$.

Proof. See Appendix A.

The above result helps us to provide $N$ equations that will be helpful for determining the vectors $\tilde{k}_{1}$, $\tilde{k}_{2}$:

Theorem 4

Suppose that the $N$ solutions $s_{i}$, $i=1,2,\ldots,N,$ with $Re(s_{i})>0$ to equation (31) are distinct. Let $\tilde{a}_{i}=(a_{i,1},\ldots,a_{i,N})^{T}$ be a nonzero column vector satisfying

[TABLE]

where $\tilde{0}$ is a row vector of zeros. Then, for $i=1,\ldots,N$, $\tilde{v}(s_{i})\tilde{a}_{i}=0$, or equivalently,

[TABLE]

where $\tilde{v}(s)=(v_{1}(s),\ldots,v_{N}(s))$ with

[TABLE]

Proof. Since $s_{i}$ satisfies equation (31), then there is a nonzero column vector $\tilde{a}_{i}$ such that

[TABLE]

Post-multiplying (28) for $s=s_{i}$ with $\tilde{a}_{i}$, we get (32), which provides $N$ linear equations that relate the elements of $\tilde{k}_{1}$, $\tilde{k}_{2}$. 2. 2.

Substituting $s=\mu_{l}$, for each $l$, $l=1,\ldots,N$ in (30), we obtain $N^{2}$ equations that relate $\Phi_{W,i}(\lambda_{j})$, $\Phi_{W,i}(\mu_{j})$, $i,j=1,\ldots,N$. More precisely, substituting $s=\mu_{l}$, $l=1,\ldots,N$ in (30) we have for $j,l=1,\ldots,N$:

[TABLE] 3. 3.

We also need $N(N-1)$ additional equations. These are obtained from (30) as follows: For each $j=1,\ldots,N$, substitute in (30) $s=\lambda_{k}$, $k\neq j$, $k=1,\ldots,N$. Then, we have for $j=1,\ldots,N$, $k\neq j$:

[TABLE]

Then, (32), (34), (35) constitute $2N^{2}$ linear equations for the $2N^{2}$ unknowns $\Phi_{W,i}(\lambda_{j})$, $\Phi_{W,i}(\mu_{j})$, $i,j=1,\ldots,N$. Once $\Phi_{W,i}(\lambda_{j})$, $\Phi_{W,i}(\mu_{j})$, $i,j=1,\ldots,N$ are known, the vectors $\tilde{k}_{1}$, $\tilde{k}_{2}$ are known, and thus, $\tilde{v}(s)$ is known, which in turn implies that $\tilde{\Phi}^{T}_{W}(s)$ is known.

4.3 Steady-state moments

Once $\tilde{k}_{1}$, $\tilde{k}_{2}$ are known, the entire transform vector $\tilde{\Phi}^{T}_{W}(s)$ is known, and ready to be used to compute the moments of the steady-state waiting time distribution. They are given in the following theorem. Denote first by:

[TABLE]

Note that $\Gamma_{0}=I$ and $\Gamma_{1}=\Gamma$. Then, we have the following result.

Theorem 5

The moment-vectors $\tilde{m}_{r}$ satisfy the following recursive equations: $\tilde{m}_{0}=\tilde{\pi}$,

[TABLE]

while for $r\geq 3$:

[TABLE]

Proof. We have that $\tilde{\Phi}^{T}_{W}(s)=\sum_{r=0}^{\infty}(-1)^{r}\tilde{m}_{r}\frac{s^{r}}{r!}$, $B^{*}(s)=\sum_{r=0}^{\infty}(-1)^{r}\Gamma_{r}\frac{s^{r}}{r!}$, thus,

[TABLE]

The equation in (28) is now written as

[TABLE]

or equivalently

[TABLE]

Equating the coefficients of $s^{0}$ we get $\tilde{m}_{0}(pP-I)=-q\tilde{\pi}$, which is true since $\tilde{m}_{0}=\tilde{\pi}$. Next, equating the coefficients of $s^{1}$, we get the expression for $\tilde{m}_{1}$ (note that since $p\in(0,1)$, $pP-I$ is invertible). We get that by equating the coefficients of $s^{2}$, $\tilde{m}_{2}$ is given in terms of $\tilde{m}_{1}$, $\tilde{m}_{0}=\tilde{\pi}$, as given in (36). Proceeding similarly, equating coefficients of $s^{r}$, $r\geq 3$ we get (37).

4.4 Asymptotic behaviour

Let $G(s):=H(s)+sI-\Lambda$, and $Adj(G(s))$ denotes its adjoint matrix. From Theorem 2 we know that

[TABLE]

where $G(s)$, is a vector of analytic functions for $Re(s)>0$. Assume that the LST $\beta_{i}^{*}(s)$ exists in a neighborhood of the origin. Then, due to (38), the functions $\Phi_{W,i}(s)$ are analytic for all $s$ with $Re(s)>-R$, where $-R$ denotes the zero of $det(G(s))$ with the largest real part in the negative halfplane. From the damping property of the Laplace transforms, we have $\mathcal{L}(e^{Rt}F_{W,i}(t))=\Phi_{W,i}(s-R)$, $i=1,\ldots,N,$ where $F_{W,i}(t)$ denotes the distribution of $W$ when the background Markov chain is in state $i$, and $\mathcal{L}(g(t))$ denotes the LST of the function $g(t)$. Therefore, given that the limit exists,

[TABLE]

where $\tilde{F}^{T}_{W}(t)=(F_{w,1}(t),\ldots,F_{W,N}(t))$. For convenience, let $s=-R$ be a simple pole of $\Phi_{W,i}(s)$, then we obtain, using de L’Hopital rule:

[TABLE]

Therefore, the waiting time distribution decays exponentially at rate $R$, and the corresponding constants are given in (39).

5 A simple numerical example

We now numerically illustrate the theoretical findings in Sections 3, 4. Let $N=2$, and $\Lambda=diag(\lambda_{1},\lambda_{2})=diag(2,3)$, $M=diag(\mu_{1}=10/u,\mu_{2}=8/u)$, $A|Z=j\sim exp(\lambda_{j})$.

Regarding Model I

Here we assume that $a=0.2$, $S|Z=j\sim exp(\mu_{j})$.

Note that in order to obtain numerical results we need to first obtain the unknowns $c_{\omega,w}$, $\omega=1,2,\ldots,l_{j}+\sum_{k=1}^{N}m_{k}$, $j=1,\ldots,N$, as stated in Theorem 1. In doing this, we firstly have to derive a system of $\sum_{j=1}^{N}l_{j}+N\sum_{k=1}^{N}m_{k}$ equations. Crucial ingredients of that system are the values of (16) for certain values of $s$; i.e., see (18), (19). Thus, a critical issue refers to the number of iterations that we need. Numerical results show that we need a (relatively) small finite number of product form terms to derive the Laplace-Stieltjes transform vector $\tilde{\Phi}_{W}(s)$. In doing that, we truncate the infinite sums of products at a specific value of $k$, such that the absolute maximum norm of the difference of two consecutive terms to be smaller than $10^{-7}$, e.g., for the derivation of $\tilde{\Phi}_{W}(a\delta_{r,j})$, let $u_{r,j}(k):=\prod_{d=0}^{k-1}R(a^{d+1}\delta_{r,j})\tilde{V}(a^{k+1}\delta_{r,j})$, $r=1,\ldots,l_{j}$, $j=1,\ldots,N$. Then, we find the maximum $k_{r,j}$ such that $||u_{r,j}(k_{r,j})-u_{r,j}(k_{r,j}+1)||\leq 10^{-7}$.

The parameter $u$ is used to explore the effect of increasing the expected service duration on the mean workload. We compare two cases: namely the case where interarrival and service times are autocorrelated i.e., where $P=\begin{pmatrix}0&1\\ 1&0\end{pmatrix}$, and the case where there is no autocorrelation, i.e., where $P=\begin{pmatrix}0.5&0.5\\ 0.5&0.5\end{pmatrix}$. In such a scenario, we have positive cross-correlation between the interarrival and service times. Figure 1 shows that the mean workload is higher when there is autocorrelation, and increases as $u$ increases.

Regarding Model II

For convenience, we assume that $S|Z=j\sim exp(\mu_{j})$, and $C|Z=j\sim exp(\theta_{j})$, $j=1,2,$ with $\theta_{1}=8$, $\theta_{2}=6$. $P=\begin{pmatrix}0.2&0.8\\ 0.6&0.4\end{pmatrix}$. Figure 2 shows the expected waiting time as a function of $p$, for increasing values of $u$. Note that as $u$ increases, the mean waiting time also increases, since the mean service time increases. That increase becomes more apparent as $p$ increases, that is, when our model is more likely to operate as a regular MAP/G/1 queue rather than an alternating service model.

6 Concluding remarks and suggestions for future research

In this work, we focused on the stationary analysis of a Markov-modulated version of a multiplicative Lindley-type recursion of the form $W_{n+1}=[V_{n}W_{n}+Y_{n}(V_{n})]^{+}$. We explicitly analyzed two cases: a) Model I: $V$ equals either 1 with probability $p_{1}$, or equals $a\in(0,1)$ with probability $p_{2}$, or it is negative with probability $p_{3}=1-p_{1}-p_{2}$, and $Y_{n}(V_{n})=S_{n}-A_{n+1}$, where $A_{n}$, $S_{n}$ have a rational LST. b) Model II: $V$ equals either one with probability $p\in(0,1)$ or minus one with probability $q=1-p$, and $Y_{n}(V_{n})$ have a general form and are also written as a difference of two nonnegative random variables that are dependent both on $V$. In both cases, $Y_{n}(V_{n})$ depend on the state of an irreducible discrete-time Markov chain with finite state space.

The analysis of Model I resulted in a vector-valued functional equation for the LST of the stationary distribution of $W$ that was solved recursively, by using the Wiener-Hopf theory. In Model II, we come up with a semi-Markov queueing system where the interarrival and service times depend linearly and randomly on the waiting times. We provided explicit stationary characteristics.

Cases which might allow explicit analysis are, for example: a) Allowing $D_{n}$ to follow a mixed Erlang distribution will not complicate the analysis in Model II. b) A challenging task might be to consider an additional dependence structure among $A_{n}$, $S_{n}$, and $D_{n}$, $C_{n}$, based on a copula. c) The analysis in Model II might be extended to the case where $V$ takes any finite number of negative values. Clearly, additional numerical experiments should be added in order to provide more insight into the stationary characteristics.

Acknowledgements

The author gratefully acknowledges the Empirikion Foundation, Athens, Greece (www.empirikion.gr) for the financial support of this work. This work is dedicated to the memory of Antonis Kontis.

Appendix A Proof of Theorem 3

The proof follows the lines in [2]. We first assume that $\beta_{i}^{*}(s)$ are analytic for all $s$, $Re(s)>-\zeta$, for some $\zeta>0$. Consider the circle $C_{\epsilon}$ with center $\widehat{\lambda}=max\{\lambda_{1},\ldots,\lambda_{N}\}$ and radius $\epsilon+\widehat{\lambda}$, $0<\epsilon<\zeta$. Let $G_{\epsilon}(s,u):=uH(s)+sI-\Lambda$, $u\in[0,1]$. We first prove that $det(G_{\epsilon}(s,u))\neq 0$, $u\in[0,1]$, $s\in C_{\epsilon}$. Indeed, $G_{\epsilon}(s,u)$, $s\in C_{\epsilon}$, $Re(s)\geq 0$ is diagonally dominant for $u\in[0,1]$ since

[TABLE]

Hence, its determinant is nonzero for $s\in C_{\epsilon}$, $Re(s)\geq 0$.

Note that $det(H(s)+sI-\lambda)=0$, is equivalent to $det(pB^{*}(s)P+s\Lambda^{-1}-I)=0$. To prove this for $s\in C_{\epsilon}$, $Re(s)<0$, first note that the determinant is nonzero if and only if $s=0$ is not an eigenvalue of $upB^{*}(s)P+s\Lambda^{-1}-I$. So in a neighbor of $s=0$ we write

[TABLE]

and we see that for $(s,u)$ close to $(0,1)$, the above matrix is a perturbation of $pP-I$ that do not has an eigenvalue 0 for $p\in(0,1)$. Thus, the determinant is nonzero for $s\in C_{\epsilon}$, $Re(s)<0$.

Let $f(u)$ be the number of zeros of $det(upB^{*}(s)P+s\Lambda^{-1}-I)$ inside $C_{\epsilon}$:

[TABLE]

Clearly $f(0)=N$, since $det(s\Lambda^{-1}-I)=\prod_{j=1}^{N}(\frac{s}{\lambda_{i}}-1)$, and thus, the zeros are $s_{i}=\lambda_{i}$ and exactly $N$ have $Re(s)\geq 0$, and thus, are inside $C_{\epsilon}$. Since $f(u)$ is an integer-valued continuous function on $[0,1]$, it is constant, therefore, $f(1)=N$. As $\epsilon\to 0$, we conclude that $det(pB^{*}(s)P\Lambda-sI-\Lambda)$ has exactly $N$ zeros inside or on $C_{0}$.