A Wavelet-Based Computational Framework for a Block-Structured Markov Chain with a Continuous Phase Variable

Abstract

We consider the computing issues of the steady probabilities for block-structured discrete-time Markov chains that are of upper-Hessenberg or lower-Hessenberg transition kernels with a continuous phase set. An effective computational framework is proposed based on the wavelet transform, which extends and modifies the arguments in the literature for quasi-birth-death (QBD) processes. A numerical procedure is developed for computing the steady probabilities based on the fast discrete wavelet transform, and several examples are presented to illustrate its effectiveness.

1. Introduction

Consider a two-dimensional block-structured discrete-time Markov chain (DTMC) $\{(L_n,X_n):n\in \mathbb{N}\}$ on the state space $\mathbb{N}\times \mathbb{R}$, where $\mathbb{N}$ and $\mathbb{R}$ are sets of non-negative integers and real numbers, respectively. Denote by $\mathcal{B}\left(\mathbb{R}\right)$ the Borel $\sigma$-field of the set $\mathbb{R}$. The transition probability law is time homogeneous and is characterized by the following transition kernel

$\begin{array}{c}\hfill P_{ij}(x,A)=\mathbb{P}\{(L_{n+1},X_{n+1})\in j\times A|(L_n,X_n)=(i,x)\},\end{array}$

where $i,j\in \mathbb{N}$, $x\in \mathbb{R}$ and $A\in \mathcal{B}\left(R\right)$. Recall that a two-dimensional function $F(x,A)$ is called a kernel if it is a measurable function in x for each $A\in \mathcal{B}\left(R\right)$, and a non-negative measure on $\mathbb{R}$ for each $x\in \mathbb{R}$. When $A=(-\infty ,y]$, we write $F(x,A)$ to be $F(x,y)$ for simplicity. Note that the kernel function $P_{ij}(x,y)$ is stochastic in the sense that $P_{ij}(x,\infty ):={lim}_{y\to \infty }{\sum }_{j\ge 0}{P}_{ij}(x,y)=1$ for all i and all x. The level and phase of each state $(i,x)$ are respectively represented by the first component i and the second component x. For any $i\ge 0$, define ${\ell }_i=\{(i,x):x\in \mathbb{R}\}$ to be the i level set. Then, the state space E can be decomposed as $E={\bigcup }_{i=0}^{\infty }{\ell }_i$. For $n\ge 1$, the corresponding n-step transition kernel is given by

$P^n(i,x;j,A)=\sum _{k\in \mathbb{N}}{\int }_{\mathbb{R}}{P}^{n-1}(i,x;k,dz)P(k,z;j,A)=\mathbb{P}\{X_n\in j\times A|X_0=(i,x)\}.$

Two different types of block-structured discrete-time Markov chains are the focus of this paper. The first one is the discrete-time GI/M/1-type Markov chain, whose transition kernel matrix $P_{GI}(x,y):={\left(P_{GI}(i,x;j,y)\right)}_{i,j\in \mathbb{N}}$ is level independent and has the following lower-Hessenberg block form:

(1)$\begin{array}{c}{P}_{GI}(x,y)=\left(\begin{array}{ccccc}{B}_0(x,y)& A_0(x,y)& 0& 0& \cdots \\ B_1(x,y)& A_1(x,y)& A_0(x,y)& 0& \cdots \\ B_2(x,y)& A_2(x,y)& A_1(x,y)& A_0(x,y)& \cdots \\ B_3(x,y)& A_3(x,y)& A_2(x,y)& A_1(x,y)& \cdots \\ ⋮& ⋮& ⋮& \ddots & \ddots \end{array}\right).\hfill \end{array}$

The second one is the discrete-time M/G/1-type Markov chain, whose transition kernel matrix $P_M(x,y):={\left(P_M(i,x;j,y)\right)}_{i,j\in \mathbb{N}}$ is level independent and has the following upper-Hessenberg block form:

(2)$P_M(x,y)=\left(\begin{array}{ccccc}{B}_0(x,y)& B_1(x,y)& B_2(x,y)& B_3(x,y)& \cdots \\ A_0(x,y)& A_1(x,y)& A_2(x,y)& A_3(x,y)& \cdots \\ 0& A_0(x,y)& A_1(x,y)& A_2(x,y)& \cdots \\ 0& 0& A_0(x,y)& A_1(x,y)& \cdots \\ ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right),$

These block-structured Markov chains are of the special features that the transition of the level is skip-free to the right or skip-free to the left, respectively.

Tweedie [1] proposed the GI/M/1-type Markov chain with a continuous phase set and demonstrated that the positive recurrent GI/M/1-type Markov chain is of the operator-geometric stationary distribution. Thus, Tweedie [1] extended the well-known results for the GI/M/1-type Markov chain with a finite phase set, which was derived by Neuts [2]. Tweedie’s finding was later applied by Breuer [3] to investigate the stationary distribution for the embedded GI/G/k queue with a Lebsegue-dominated inter-arrival time distribution. A positive recurrent tridiagonal block-structured quasi-birth-death (QBD) process with a continuous phase set, as well as a computational framework of its stationary distribution, are investigated by Nielsen and Ramaswami [4]. They also demonstrated the motivation for investigating a model with a continuous phase set. The computational framework was recently extended and improved by Jiang et al. [5] by incorporating the wavelet transform approach.

The GI/M/1-type and M/G/1-type Markov chains with a finite phase set were investigated systematically by Neuts in 1981 [2] and 1989 [6], respectively. Effective solver tools for solving the stationary distribution for these chains were developed by Bini et al. in [7], based on the algorithms collected in [8]. It is known that the matrices R and G are key matrices for solving stationary distributions for GI/M/1-type and M/G/1-type Markov chains, respectively. Since R and G are closely connected by Ramaswami dual and Bright dual, the computation of matrix R for GI/M/1-type chains can be reduced to the computation of matrix G for M/G/1-type Markov chains ([9,10,11,12]). Several effective algorithms have been developed to compute the matrix G, such as functional iteration, Newton iteration, invariant subspace method, cyclic reduction and Ramaswami Reduction. Please refer to [13] for a detailed description of the algorithms.

As far as we know, the following two issues are still not well addressed in the literature:

(i) For a positive recurrent GI/M/1-type Markov chain with a continuous phase set, numerical algorithms for computing the stationary distribution are missing, although the theoretical framework has been established in [1],

(ii) M/G/1-type Markov chains are of the same importance as GI/M/1-type Markov chains. However, both the theoretical and computational framework are missing for M/G/1-type Markov chains with a continuous phase set.

The current research is motivated to investigate the above two issues. This paper is organized into six sections. We provide an overview of DTMCs on a general state space and the wavelet series expansion in two dimensions in Section 2. The GI/M/1-type Markov chains are introduced in Section 3, most of which are well known in the literature [1], except for the computational analysis. The analysis of stationary distributions for M/G/1-type Markov chains is performed in Section 4. Numerical experiments, including a brief description of numerical algorithms and two illustrative examples, are presented in Section 5. Comparisons among different algorithms are executed with respect to the accuracy and speed of calculation. Conclusions are presented in Section 6. Please refer to Table A1 for a summary of frequently used notations.

2. Preliminaries

2.1. Basics about DTMCs on a General State Space

We present some basic concepts for DTMCs on a general state space. Please refer to [14] for more details.

Let ${\mathrm{\Phi }}_n$ be a DTMC on a general state space E endowed with the countably generated $\sigma$-field $\mathcal{B}\left(E\right)$. Define ${\tau }_A=\{n\ge 1:{\mathrm{\Phi }}_n\in A\}$ to be the first return time on A. For a non-negative nontrivial measure $\psi$, the chain ${\mathrm{\Phi }}_n$ is called $\psi$-irreducible if there exits a non-negative nontrivial measure $\phi$, such that ${\mathrm{\Phi }}_n$ is $\phi$-irreducible, i.e.,

$L(x,A):=P\{{\tau }_A<\infty |{\mathrm{\Phi }}_0=x\}>0$

for any $A\in \mathcal{B}\left(E\right)$, $\phi \left(A\right)>0$ and any $x\in A$, and $\psi$ is a maximal irreducible measure with respect to $\phi$. A set $A\in \mathcal{B}\left(E\right)$ is called a Harris recurrent if $L(x,A)=1$ for all $x\in A$. The chain ${\mathrm{\Phi }}_n$ is called a Harris recurrent if it is $\psi$-irreducible and every set in ${\mathcal{B}}_{+}\left(E\right):=\{A\in \mathcal{B}\left(E\right):\psi \left(A\right)>0\}$ is Harris recurrent. A Harris recurrent chain has an unique invariant measure $\mathrm{\Pi }$ such that

$\mathrm{\Pi }\left(A\right)={\int }_E\mathrm{\Pi }\left(dx\right)P(x,A).$

A Harris recurrent chain with a finite $\mathrm{\Pi }\left(E\right)$ is said to be Harris positive recurrent. If ${\mathrm{\Phi }}_n$ is Harris positive recurrent and aperiodic, then

$\mathrm{\Pi }\left(A\right)=\underset{n\to \infty }{lim}{P}^n(x,A),$

which implies that the limit of the transition kernel exists independently of the initial state $(i,x)$. In this case, $\mathrm{\Pi }\left(A\right)$ is called the invariant probability measure or the stationary probability distribution.

We now introduce the censored Markov chain, which will be used later to deal with the invariant probability distributions for block-structured Harris positive recurrent chains. Let A be a non-empty subset in $\mathcal{B}\left(E\right)$. Let ${\theta }_k$ be the $kth$ time that ${\mathrm{\Phi }}_n$ successively visits a state in A, i.e., ${\theta }_0=inf\left\{m\ge 0:{\mathrm{\Phi }}_m\in A\right\}$ and ${\theta }_{k+1}=inf\left\{m\ge {\theta }_k+1:{\mathrm{\Phi }}_m\in A\right\}$. The censored Markov chain ${\mathrm{\Phi }}^A=\left\{{\mathrm{\Phi }}_k^A,k\ge 0\right\}$ on A is defined by ${\mathrm{\Phi }}_k^A={\mathrm{\Phi }}_{{\theta }_k},k\ge 0$, whose one-step transition kernel is denoted by $P^A(x,B)$, $x\in E,B\in \mathcal{B}\left(E\right)$. Define

${}_A{P}^n(x,B)=P\{{\mathrm{\Phi }}_n\in B,{\mathrm{\Phi }}_m\notin A,1\le m\le n\mid {\mathrm{\Phi }}_0=x\},$

and

$U_A(x,B)=\sum _{n=1}^{\infty }{}_A{P}^n(x,B).$

When starting with ${\mathrm{\Phi }}_0=x\in A$ and $B\subseteq A$, the censored chain ${\mathrm{\Phi }}^A$ evolves according to the transition law

$P^A(x,B)=U_A(x,B)=P\{{\mathrm{\Phi }}_{{\tau }_A}\in B\mid {\mathrm{\Phi }}_0=x\}.$

2.2. Basics about Wavelet in Two Dimensions

This section is concerned on some basics about the wavelet, most of which is taken from [5] directly. Please refer to [15,16] for more details about the wavelet analysis.

Respectively denote the scaling function and the wavelet function by $\varphi$ and $\psi$. For all $j\in \mathbb{Z}$, let ${\varphi }_{j,n}\left(x\right)=2^{-j/2}\varphi \left(2^{-j}(x-2^{j}n)\right)$ and ${\psi }_{j,n}\left(x\right)=2^{-j/2}\psi \left(2^{-j}(x-2^{j}n)\right)$. Define three wavelets

$W^{\left(1\right)}(x_1,x_2)=\varphi \left(x_1\right)\psi \left(x_2\right),W^{\left(2\right)}(x_1,x_2)=\psi \left(x_1\right)\varphi \left(x_2\right),W^{\left(3\right)}(x_1,x_2)=\psi \left(x_1\right)\psi \left(x_2\right),$

and for all j in $\mathbb{Z}$,

$W_{j,n_1,n_2}^{\left(k\right)}(x_1,x_2)=\frac{1}{2^j}{W}^{\left(k\right)}\left(\frac{x_1-2^j{n}_1}{2^j},\frac{x_2-2^j{n}_2}{2^j}\right),n_1,n_2\in \mathbb{Z},1\le k\le 3.$

Now, we consider the wavelet series expansion of a two-dimensional function. For each $i\in \mathbb{Z}$, define column vectors ${\varphi }_i\left(x\right)=[{\varphi }_{i,n}\left(x\right):n\in \mathbb{Z}]$, ${\psi }_i\left(x\right)=[{\psi }_{i,n}\left(x\right):n\in \mathbb{Z}]$, and ${\zeta }_i\left(x\right)={\left[{\varphi }_i^T\left(x\right),{\psi }_i^T\left(x\right)\right]}^T$. By Lemma 3.1 in [5], any function $u(x,y)\in {\mathit{L}}^2\left({\mathbb{R}}^2\right)$ can be expanded as follows

(3)$u(x,y)={\zeta }^T\left(x\right)\overline{U}\zeta \left(y\right),$

where $\zeta \left(x\right)=[{\zeta }_i\left(x\right):i\in \mathbb{Z}]$ is a column vector, the diagonal blocks of $\overline{U}$ are written as

$\begin{array}{c}\hfill {\overline{U}}_i=\left[\begin{array}{cc}0\hfill & {\overline{U}}_i^{\left(1\right)}\hfill \\ {\overline{U}}_i^{\left(3\right)}\hfill & {\overline{U}}_i^{\left(3\right)}\hfill \end{array}\right]\end{array}$

with ${\left({\overline{U}}_i^{\left(k\right)}\right)}_{m,n}=<u,W_{i,m,n}^{\left(k\right)}>$.

Let $U(x,y)$ be a kernel function whose density function is assumed to exist and is denoted by $u(x,y):=\frac{\partial U(x,y)}{\partial y}$. On the one hand, by performing (3) for the density of the kernel function $U(x,y)$, which is referred to as the wavelet transform (WT), we can find the matrix $\overline{U}$, also known as the associated matrix of $U(x,y)$. On the other hand, for a given associated matrix $\overline{U}$, we can find the density function $u(x,y)$ by performing (3) in the other side, which is called the inverse wavelet transform (IWT).

As you will see in the following sections, it is crucial to deal with the convolution operations of the transition kernels in order to investigate the stationary distributions. The wavelet transform is introduced to transform these convolution operations into matrix operations by expanding the kernels using the wavelet series. For any $A\in \mathcal{B}\left(\mathcal{E}\right)$, define the convolution $C_1*C_2$ of two kernel functions $C_1(x,y)$ and $C_2(x,y)$ by

$\begin{array}{c}\hfill C_1*C_2(x,A)={\int }_{\mathbb{R}}{C}_1(x,dz)C_2(z,A),\end{array}$

and define $C_1^{\left(k\right)}(x,y)$ recursively by

$\begin{array}{c}\hfill C_1^{\left(k\right)}(x,y)={\int }_{\mathbb{R}}{C}_1^{(k-1)}(x,dz)C_1(z,y)={\int }_{\mathbb{R}}{C}_1(x,dz)C_1^{(k-1)}(z,y),\end{array}$

where $C_1^{\left(0\right)}(x,x)=1$ and $C_1^{\left(0\right)}(x,y)=0$ for any $y\ne x$. If $\nu$ is a signed measure on E, we write

$\nu *C_1\left(A\right)={\int }_E\nu \left(dx\right)C_1(x,A).$

In order to expand the kernel functions through the wavelet transform, we need the following assumption and theorem, which are both taken from [5].

Assumption 1.

All the kernel functions $B_i(x,y)$ and $A_i(x,y)$, $i\ge -1$ belong to ${\sum }_H$, where $H\subset \mathbb{R}$ is of finite Lebesgue measure, and ${\sum }_H$ is the set of kernel functions $U(x,y)$ having a density function $u(x,y)$ equaling to zeros outside of $H\times H$.

Theorem 1

([5]). Let $\{F_k(x,y),k\ge 1\}$ be a sequence of kernel functions in ${\sum }_H$. Denote their density functions by $\{f_k(x,y),k\ge 1\}$, and their associated matrices by $\{{\overline{F}}_k,k\ge 1\}$.

(i) For any fixed n, the convolution kernel function $F_1*F_2*\cdots *F_n(x,y)$ is also in ${\sum }_H$, and its associated matrix is ${\prod }_{k=1}^n{\overline{F}}_k$;

(ii) For any fixed n, the additive kernel function $(F_1+F_2+\cdots +F_n)(x,y)$ is also in ${\sum }_H$ and its associated matrix is ${\sum }_{k=1}^n{\overline{F}}_k$;

(iii) If $f_n(x,y)$ converges to $f(x,y)$, then the kernel function $F(x,y):={\int }_{-\infty }^{y}f(x,y)dy$ is also in ${\sum }_H$, and its associated matrix is $\overline{F}={lim}_{n\to \infty }{\overline{F}}_n$.

3. $\mathit{GI}/\mathit{M}/\mathbf{1}$-Type Markov Chains

Consider a GI/M/1-type Markov chain $(L_n,X_n)$, whose transition law P given by (1) satisfies that for any $C\in \mathcal{B}\left(\mathbb{R}\right)$

$P(i,x;j,\mathbb{R})=0,j>i+1,$

$P(i,x;j,C)=A_{i-j+1}(x,C),i\ge j-1,j\ge 1,$

$P(i,x;0,C)=B_i(x,C),i\ge 0.$

Define the kernel $R(x,C)$ to be the expected number of visits to $(i+1)\times C$, starting from $(i,x)$ under the taboo set of ${\bigcup }_{k=0}^i{\ell }_k$. From [1], we know that the censored Markov chain ${(L_n,X_n)}^{{\ell }_0}$ of the GI/M/1-type Markov chain on the zero level set ${\ell }_0$ has the following transition kernel

$P_{GI}^{{\ell }_0}(x,C)=\sum _{k=0}^{\infty }{R}^{\left(k\right)}*B_k(x,C).$

The following theorem is taken from [1], which characterizes the invariant probability measure for GI/M/1-type Markov chains with a continuous phase set.

Theorem 2

([1]). Suppose that the GI/M/1-type Markov Chain $(L_n,X_n)$ with a continuous phase set is Ψ-irreducible and Harris positive recurrent. Then, its unique stationary probability measure Π, decomposed by $\mathrm{\Pi }\left(A\right)=({\mathrm{\Pi }}_0\left(A\right),{\mathrm{\Pi }}_1\left(A\right),{\mathrm{\Pi }}_2\left(A\right),\cdots )$, satisfies the following recursive formula

${\mathrm{\Pi }}_k\left(A\right)={\mathrm{\Pi }}_0*R^{\left(k\right)}\left(A\right),$

where the kernel $R(x,A)$ is the minimal non-negative solution of the following equation

$R(x,A)=\sum _{i=0}^{\infty }{R}^{\left(i\right)}*A_i(x,A),$

and ${\mathrm{\Pi }}_0\left(A\right)$ is uniquely determined by

${\mathrm{\Pi }}_0*P_{GI}^{{\ell }_0}(x,A)={\mathrm{\Pi }}_0\left(A\right),\sum _{k=0}^{\infty }{\mathrm{\Pi }}_0*R^{\left(k\right)}(x,E)=1.$

Applying Theorem 2 and Theorem 1, we can obtain the following theorem directly.

Theorem 3.

Suppose that the GI/M/1-type Markov Chain $(L_n,X_n)$ with a continuous phase set is ψ-irreducible and Harris positive recurrent and that Assumption 1 holds.

(i) The kernels $P_{GI}^{{\ell }_0}(x,y)$ and $R(x,y)$ are in ${\sum }_H$, whose associated matrices are, respectively, denoted by ${\overline{P}}_{GI}^{{\ell }_0}$, $\overline{R}$ and ${\overline{B}}_k$.

(ii) The invariant probability measure ${\mathrm{\Pi }}_k$ is in ${\sum }_H$. Let ${\overline{\mathrm{\Pi }}}_k$ be the associated row vector of ${\mathrm{\Pi }}_k\left(y\right)$, i.e., ${\pi }_k\left(y\right)={\overline{\mathrm{\Pi }}}_k\zeta \left(y\right)$, where ${\pi }_k\left(y\right)$ is the density of ${\mathrm{\Pi }}_k\left(y\right)$, and $\zeta \left(y\right)$ is defined in Section 2.2. Then, we have

${\overline{\mathrm{\Pi }}}_k={\overline{\mathrm{\Pi }}}_0{\overline{R}}^{\left(k\right)},$

and

${\overline{\mathrm{\Pi }}}_0{\overline{P}}_{GI}^{{\ell }_0}={\overline{\mathrm{\Pi }}}_0,{\overline{\mathrm{\Pi }}}_0{(I-\overline{R})}^{-1}\mathbf{1}=1,$

where ${\overline{P}}_{GI}^{{\ell }_0}={\sum }_{k=0}^{\infty }{\overline{R}}^{\left(k\right)}{\overline{B}}_k$ and $\mathbf{1}$ is the vector of 1’s with an appropriate dimension.

4. $\mathit{M}/\mathit{G}/\mathbf{1}$-Type Markov Chains

In this section, we consider a M/G/1-type Markov Chain $(L_n,X_n)$, whose transition law P, given by (1), satisfies that for any $C\in \mathcal{B}\left(\mathbb{R}\right)$

$P(i,x;j,C)=0,forj<i-1,i\ge 1,$

$P(0,x;j,C)=B_j(x,C),$

$P(i,x;j,C)=A_{j-i+1}(x,C),forj\ge i-1,i\ge 1.$

Define ${\tau }_{{\ell }_i}=inf\left\{n\ge 1:L_n\in {\ell }_i\right\}$ to be the first return time to the level set ${\ell }_i$ for any $i\ge 0$. For any $x\in \mathbb{R}$ and any $A\in \mathcal{B}\left(\mathbb{R}\right)$, define the following kernel function

$G(x,A)=P\left\{{\tau }_{{\ell }_i}<\infty ,X_{{\tau }_{{\ell }_i}}\in A\mid L_0=i+1,X_0=x\right\},$

which is independent of i due to the level independent structure of the chain. The first result is about the kernel $G(x,A)$, which plays a key role in analyzing M/G/1-type Markov chains.

Theorem 4.

Suppose that the M/G/1-type Markov chain $(L_n,X_n)$ is ψ-irreducible. For any $A\in \mathcal{B}\left(\mathbb{R}\right)$, the kernel $G(x,A)$ is the minimal nonnegative solution of the following equation

(4) $G(x,A)=\sum _{i=0}^{\infty }{A}_i*G^{\left(i\right)}(x,A),$

where $G^{\left(i\right)}(x,A)$ is the i-fold convolution of the kernel $G(x,A)$ itself.

Proof. 

We first show that the kernel $G(x,A)$ is a solution of Equation (4).

By conditioning on the state of the first transition, the kernel $G(x,A)$ can be decomposed as follows

(5)$\begin{array}{ccc}\hfill G(x,A)& =& {\int }_{\mathbb{R}}\sum _{i\in \mathbb{N}}\left[P\left\{L_1=i,X_1\in dy\mid L_0=1,X_0=x\right\}\right.\hfill \\ & & \times \left.P\left\{{\tau }_{{\ell }_0}<\infty ,X_{{\tau }_{{\ell }_0}}\in A\mid L_0=i,X_0=y\right\}\right]\hfill \\ & =& {\int }_{\mathbb{R}}\sum _{i\in \mathbb{N}}{A}_i(x,dy)P\left\{{\tau }_{{\ell }_0}<\infty ,X_{{\tau }_{{\ell }_0}}\in A\mid L_0=i,X_0=y\right\}.\hfill \end{array}$

We will use the inductive arguments to show

(6)$G^{\left(i\right)}(y,A)=P\left\{{\tau }_{{\ell }_0}<\infty ,X_{{\tau }_{{\ell }_0}}\in A\mid L_0=i,X_0=y\right\},i\ge 1.$

Since the chain is level independent, when $i=1$, we have

$G^{\left(1\right)}(y,A)=G(y,A)=P\left\{{\tau }_{{\ell }_k}<\infty ,X_{{\tau }_{{\ell }_k}}\in A\mid L_0=k+1,X_0=y\right\}$

for any $k\ge 0$. Suppose that $G^{\left(n\right)}$ satisfies

$G^{\left(n\right)}(y,A)=P\left\{{\tau }_{{\ell }_0}<\infty ,X_{{\tau }_{{\ell }_0}}\in A\mid L_0=n,X_0=y\right\}.$

By conditioning on the state of the first hitting on level ${\ell }_n$ and using the strong Markov property, we have

$\begin{array}{cc}& P\left\{{\tau }_{{\ell }_0}<\infty ,X_{{\tau }_{{\ell }_0}}\in A\mid L_0=n+1,X_0=y\right\}\hfill \\ & ={\int }_{\mathbb{R}}P\left\{{\tau }_{{\ell }_n}<\infty ,X_{{\tau }_{{\ell }_n}}\in dx\mid L_0=n+1,X_0=y\right\}P\left\{{\tau }_{{\ell }_0}<\infty ,X_{{\tau }_{{\ell }_0}}\in A\mid L_{{\tau }_{{\ell }_n}}=n,X_{{\tau }_{{\ell }_n}}=x\right\}\hfill \\ & ={\int }_{\mathbb{R}}P\left\{{\tau }_{{\ell }_n}<\infty ,X_{{\tau }_{{\ell }_n}}\in dx\mid L_0=n+1,X_0=y\right\}P\left\{{\tau }_{{\ell }_0}<\infty ,X_{{\tau }_{{\ell }_0}}\in A\mid L_0=n,X_0=x\right\}\hfill \\ & ={\int }_{\mathbb{R}}{G}^{\left(1\right)}(y,dx)G^{\left(n\right)}(x,A)\hfill \\ & =G^{(n+1)}(y,A).\hfill \end{array}$

Substituting (6) into (5), we have

$G(x,A)={\int }_{\mathbb{R}}\sum _{i=0}^{\infty }{A}_i(x,dy)G^{\left(i\right)}(y,A)=\sum _{i=0}^{\infty }{A}_i*G^{\left(i\right)}(y,A),$

where we exchange the order between integration and summation by Fubini theorem.

Next we demonstrate that $G(x,A)$ is the minimal non-negative solution of (4). We divide the proof into two steps.

We first define a sequence of kernels $\{T_N(x,A),N\ge 1\}$ by setting $T_0(x,A)=0$, and

$T_{N+1}(x,A)=\sum _{i=0}^{\infty }{A}_i*T_N^{\left(i\right)}(x,A),N\ge 0.$

Let $\widehat{G}(x,A)$ be any solution of Equation (4). Obviously, $\widehat{G}(x,A)\ge 0=T_0(x,A)$. Suppose that $T_{N-1}(x,A)\le \widehat{G}(x,A)$, then $T_{N-1}^{\left(i\right)}(x,A)\le {\widehat{G}}^{\left(i\right)}(x,A)$ for $i\ge 1$. Moreover, we have

(7)$T_N(x,A)=\sum _{i=0}^{\infty }{A}_i*T_{N-1}^{\left(i\right)}(x,A)\le \sum _{i=0}^{\infty }{A}_i*{\widehat{G}}^{\left(i\right)}(x,A)=\widehat{G}(x,A).$

Similarly, if we assume inductively that $T_{N-1}(x,A)\le T_N(x,A)$, we have

$T_N(x,A)=\sum _{i=0}^{\infty }{A}_i*T_{N-1}^{\left(i\right)}(x,A)\le \sum _{i=0}^{\infty }{A}_i*T_N^{\left(i\right)}(x,A)=T_{N+1}(x,A),$

and so $T_N(x,A)$ is monotonically increasing in N. Hence, the limit $T_{*}(x,A):={lim}_{N\to \infty }\uparrow T_N(x,A)$ exists. Further, we have

$T_N^{\left(k\right)}(x,A)\uparrow T_{*}^{\left(k\right)}(x,A),k\ge 1,$

By taking the limit of both sides of Equation (6) and using the dominated convergence theorem, we know that the kernel $T_{*}(x,A)$ is a solution of (4), i.e.,

$T_{*}(x,A)=\sum _{i=0}^{\infty }{A}_i*T_{*}^{\left(i\right)}(x,A),$

We further have that $T_{*}(x,A)$ is the minimal solution since $T_{*}(x,A)\le \widehat{G}(x,A)$.

Next, we need to prove that $T_{*}(x,A)=G(x,A)$. Define

$G_N(x,A)=P\left\{{\tau }_{{\ell }_0}\le N,X_{{\tau }_{{\ell }_0}}\in A\mid L_0=1,X_0=x\right\},N\ge 1.$

Obviously, we know that $G_N(x,A)\uparrow G(x,A)$. By conditioning on the state of the first transition, we have

(8)$\begin{array}{ccc}\hfill G_{N+1}(x,A)& =& {\int }_{\mathbb{R}}\sum _{i=0}^{N}P\left\{{\tau }_{{\ell }_0}\le N,X_{{\tau }_{{\ell }_0}}\in A\mid L_0=i,X_0=y\right\}\hfill \\ & \times & P\left\{L_1=i,X_1\in dy\mid L_0=1,X_0=x\right\}.\hfill \end{array}$

Denote

$M_N^{\left(i\right)}(y,A)=P\left\{{\tau }_{{\ell }_0}\le N,X_{{\tau }_{{\ell }_0}}\in A\mid L_0=i,X_0=y\right\}.$

By conditioning on the state of the first return time to level ${\ell }_{i-1}$ and repeating the same arguments, we have

(9)$\begin{array}{ccc}\hfill M_N^{\left(i\right)}(y,A)& \le & {\int }_{\mathbb{R}}P\left\{{\tau }_{{\ell }_{i-1}}\le N,X_{{\tau }_{{\ell }_{i-1}}}\in dx\mid L_0=i,X_0=y\right\}{M}_N^{(i-1)}(x,A)\hfill \\ & =& G_N^{\left(1\right)}*M_N^{(i-1)}(y,A)\hfill \\ & \le & G_N^{\left(2\right)}*M_N^{(i-2)}(y,A)\hfill \\ & \le & \cdots \hfill \\ & \le & G_N^{(i-1)}*M_N^{\left(1\right)}(y,A)\hfill \\ & =& G_N^{\left(i\right)}(y,A).\hfill \end{array}$

By (8) and (9), we can deduce that

$\begin{array}{ccc}\hfill G_{N+1}(x,A)& \le & {\int }_{\mathbb{R}}\sum _{i=0}^N{A}_i(x,dy)G_N^{\left(i\right)}(y,A)\hfill \\ & \le & {\int }_{\mathbb{R}}\sum _{i=0}^{\infty }{A}_i(x,dy)G_N^{\left(i\right)}(y,A)\hfill \\ & =& \sum _{i=0}^{\infty }{A}_i*G_N^{\left(i\right)}(x,A)\hfill \end{array}$

Finally, note that $G_1(x,A)=A_0(x,A)=T_1(x,A)$, and so from (6), we have by induction $G_N(x,A)\le T_N(x,A)$. Taking the limit as $N\to \infty$ gives $G(x,A)\le T_{*}(x,A)$, as required. □

In the following, we will investigate numerical computing issues of the invariant probability distribution for M/G/1-type chains. The key point is to set up the Ramaswami algorithm, a well-known result for M/G/1-type chains with finite phases, for the M/G/1-type chains with continuous phases.

Theorem 5.

Suppose that the M/G/1-type Markov chain $(L_n,X_n)$ is ψ-irreducible and Harris positive recurrent. Let the unique invariant probability measure be Π with $\mathrm{\Pi }\left(C\right)=({\mathrm{\Pi }}_0\left(C\right),{\mathrm{\Pi }}_1\left(C\right),$${\mathrm{\Pi }}_2\left(C\right),\cdots )$, $C\in \mathcal{B}\left(\mathbb{R}\right)$. Then, the measure Π satisfies the following recursive formula

(10)${\mathrm{\Pi }}_k\left(C\right)={\mathrm{\Pi }}_0*{\widehat{B}}_k\left(C\right)+\sum _{i=1}^k{\mathrm{\Pi }}_i*{\widehat{A}}_{k+1-i}\left(C\right),$

where

${\widehat{B}}_k(x,C)=\sum _{i=k}^{\infty }{B}_i*G^{(i-k)}(x,C),k\in \mathbb{N},$

(11)${\widehat{A}}_m(x,C)=\sum _{i=m}^{\infty }{A}_i*G^{(i-m)}(x,C),1\le m\le k,$

and ${\mathrm{\Pi }}_0$ is a unique solution of the equation ${\mathrm{\Pi }}_0\left(C\right)={\mathrm{\Pi }}_0*{\widehat{B}}_0\left(C\right)$.

Proof. 

By (6), we know that for any $v\ge 1$ and $i\ge 0$, the kernel function $G^{\left(v\right)}(x,C)$ is the probability that the Markov chain first returns to level ${\ell }_i$ by hitting the state $(i,C)$, given that it starts from the state $(i+v,x)$. The transition kernel function of the Markov chain embedded at epochs of visits to the set $A={\bigcup }_{m=0}^{\infty }{\ell }_m$ is given by

$P^A(x,C)=\left(\begin{array}{ccccc}{B}_0(x,C)& B_1(x,C)& \cdots & B_{k-1}(x,C)& {\widehat{B}}_k(x,C)\\ A_0(x,C)& A_1(x,C)& \cdots & A_{k-1}(x,C)& {\widehat{A}}_k(x,C)\\ 0& A_0(x,C)& \cdots & A_{k-2}(x,C)& {\widehat{A}}_{k-1}(x,C)\\ ⋮& ⋮& & ⋮& ⋮\\ 0& 0& \cdots & A_0(x,C)& {\widehat{A}}_1(x,C)\end{array}\right).$

We now explain how to determine the transition kernel $P^A(x,C)$. The first k block columns of the kernel function $P^A(x,C)$ are the same as those of $P(x,C)$, since the chain $(L_n,X_n)$ can only move down by one level at a time. As for the $(k+1)$th (i.e., last) block column of $P^A(x,C)$, its first entry is as follows.

$\begin{array}{ccc}\hfill {\widehat{B}}_k(x,C)& =& P\{{(L_1,X_1)}^A\in k\times C\mid {(L_0,X_0)}^A=(0,x)\}\hfill \\ & =& P\{L_1=k,X_1\in C\mid L_0=0,X_0=x\}\hfill \\ & & +\sum _{i=k+1}^{\infty }{\int }_{\mathbb{R}}[P\{L_1=i,X_1\in dy\mid L_0=0,X_0=x\}\hfill \\ & & \times P\{{\tau }_{{\ell }_k}<\infty ,X_{{\tau }_{{\ell }_k}}\in C\mid L_1=i,X_1=y\}]\hfill \\ & =& B_k(x,C)+\sum _{i=k+1}^{\infty }{B}_i*G^{(i-k)}(x,C)\hfill \\ & =& \sum _{i=k}^{\infty }{B}_i*G^{(i-k)}(x,C).\hfill \end{array}$

The equality (11) can be proved in a similar way.

Since this chain $(L_n,X_n)$ is $\psi$-irreducible and a Harris positive recurrent, for $\forall x\in X,M\subseteq E$, starting from x, the set M will almost certainly be returned infinitely, and so is the censored Markov chain ${(L_n,X_n)}^A$. Thus, ${(L_n,X_n)}^A$ is also $\psi$-irreducible and a Harris positive recurrent. Let ${\mathrm{\Pi }}^A\left(C\right)=({\mathrm{\Pi }}_0^A\left(C\right),{\mathrm{\Pi }}_1^A\left(C\right),\cdots ,{\mathrm{\Pi }}_k^A\left(C\right))$ be the unique invariant probability measure of ${(L_n,X_n)}^A$.

Next, we will demonstrate that $({\mathrm{\Pi }}_0\left(C\right),{\mathrm{\Pi }}_1\left(C\right),\cdots ,{\mathrm{\Pi }}_k\left(C\right))$ is also an invariant measure of the censored chain ${\mathrm{\Phi }}_n^A$. Define the measure ${\mathrm{\Pi }}^{\circ }$ by

${\mathrm{\Pi }}_i^{\circ }\left(C\right):={\int }_{\mathbb{R}}{\mathrm{\Pi }}_i^A\left(dx\right)U_A(x,i\times C),i\in N.$

By Propostion 10.4.8 in [14], we know that

(12)${\mathrm{\Pi }}_i^A\left(C\right)={\mathrm{\Pi }}_i^{\circ }\left(C\right),0\le i\le k$

and that ${\mathrm{\Pi }}^{\circ }$ is invariant measure for $(L_n,X_n)$. Since $(L_n,X_n)$ is assumed to be a Harris positive recurrent, the invariant measure is unique up to a constant. This shows that ${\mathrm{\Pi }}^{\circ }\left(C\right)=c\mathrm{\Pi }\left(C\right)$ for some constant c, from which and (12), and we have

${\mathrm{\Pi }}_i^A\left(C\right)=c{\mathrm{\Pi }}_i\left(C\right),\mathrm{for}0\le i\le k.$

Since ${\sum }_{i=0}^k{\mathrm{\Pi }}_i^A(-\infty ,\infty )=1$, we can obtain

$c=\frac{1}{{\sum }_{i=0}^k{\mathrm{\Pi }}_i(-\infty ,\infty )}.$

Thus, we have proved that $({\mathrm{\Pi }}_0,{\mathrm{\Pi }}_1,{\mathrm{\Pi }}_2,\cdots ,{\mathrm{\Pi }}_k)$ is an invariant measure of ${\mathrm{\Phi }}^A$. Taking into account the last block equation of ${\mathrm{\Pi }}^A\left(C\right)=\mathrm{\Pi }*P^A\left(C\right)$, we have

(13)${\mathrm{\Pi }}_k\left(y\right)={\mathrm{\Pi }}_0*{\widehat{B}}_k\left(y\right)+\sum _{i=1}^k{\mathrm{\Pi }}_i*{\widehat{A}}_{k+1-i}\left(y\right),k\ge 1.$

This proves (10).

To determine ${\mathrm{\Pi }}_0\left(C\right)$, we reset $A={\ell }_0$ and consider the censored chain ${(L_n,X_n)}^A$, whose transition kernel is given by

$P^{{\ell }_0}(x,C)={\widehat{B}}_0(x,C).$

By (13), we know that ${\mathrm{\Pi }}_0\left(C\right)={\mathrm{\Pi }}_0*{\widehat{B}}_0\left(C\right)$. □

Applying Theorem 5 and performing the wavelet series expansion, we can obtain the following theorem directly.

Theorem 6.

Suppose that the M/G/1-type Markov chain $(L_n,X_n)$ is ψ-irreducible and Harris positive recurrent and that Assumption 1 holds.

(i) The kernels $G(x,y)$, ${\overline{A}}_k(x,y)$ and ${\overline{B}}_k(x,y)$ are in ${\sum }_H$, whose associated matrices satisfies that

${\overline{\widehat{B}}}_k=\sum _{i=k}^{\infty }{B}_i{\overline{G}}^{i-k},{\overline{\widehat{A}}}_k=\sum _{i=k}^{\infty }{A}_i{\overline{G}}^{i-k}.$

(ii) The invariant probability measure ${\mathrm{\Pi }}_k$ is in ${\sum }_H$. Let ${\overline{\mathrm{\Pi }}}_k$ be the associated row vector of ${\mathrm{\Pi }}_k\left(y\right)$. Then, the associated matrices satisfy

(14) ${\overline{\mathrm{\Pi }}}_k=\left[{\overline{\mathrm{\Pi }}}_0{\overline{\widehat{B}}}_k+\sum _{i=1}^{k-1}{\overline{\mathrm{\Pi }}}_i{\overline{\widehat{A}}}_{k+1-i}\right]{\left(I-{\overline{\widehat{A}}}_1\right)}^{-1},k\ge 1,$

where ${\overline{\mathrm{\Pi }}}_0={\overline{\mathrm{\Pi }}}_0{\overline{\widehat{B}}}_0$.

Remark 1.

(i) We note that the entries in the associated matrices of a kernel function may be negative. Hence, the associated transition kernel matrices ${\overline{A}}_k^{\prime }s$ and ${\overline{B}}_k^{\prime }s$ cannot construct a stochastic transition matrix.

(ii) We now consider numerical algorithms for computing the associated matrix $\overline{G}$. In the literature, several known algorithms, including the functional iteration ([17,18]), Newton iteration ([19]), invariant subspace method ([20]), cyclic reduction ([21]) and Ramaswami Reduction ([22]), have been developed to solve the G-matrix for M/G/1-type chains with a finite phase. For a collection of these algorithms, please refer to [13]. Similar to what we did in Theorems 4 and 5, we can set up the corresponding algorithms for $\overline{G}$ by modifying these algorithms from the finite phase to the general phase. We omit the details in order to avoid tedious presentations.

5. Numerical Experiments

5.1. Discrete Wavelet Transforms

We need to perform discrete wavelet transforms for numerical experiments. Without a loss of generality, we assume that the phase space is taken to be $\mathbb{R}$. In the following, we only give a simple presentation of the computation framework; please refer to Section 5 in [5] for more details.

We first consider numerical issues of the M/G/1-type Markov chain with a continuous phase, which are divided into the following steps:

Step 1: Choose appropriate real numbers $\underline{y}$, $\overline{y}$ and positive integer N. Then, evenly sample N points from the truncated phase space $[\underline{y},\overline{y}]$. Performing the DWT in Algorithm 5.1 in [5] to kernels $A_i(x,y)$ and $B_i(x,y)$ produces the associated sample matrices ${\left(A_i\right)}_{asm}$ and ${\left(B_i\right)}_{asm}$.

Step 2: Solve the the associated sample matrix $G_{asm}$ through the algorithms listed in (ii) of Remark 1, such as functional iteration, Newton iteration, invariant subspace method, cyclic reduction and Ramaswami Reduction.

Step 3: Solve the associated sample invariant probability vector ${\mathrm{\Pi }}_{asm}$ using Theorem 6.

Step 4: Performing the IDWT in Algorithm 5.2 in [5] to the matrix of $G_{asm}$ and the vector ${\mathrm{\Pi }}_{asm}$ produces the kernels $G(x,y)$ and $\mathrm{\Pi }\left(x\right)$.

Now, we consider numerical issues of GI/M/1-type Markov chains with the continuous phase, which are also divided into four steps. The first step and the last step are the same as that for the M/G/1-type Markov chains. In Step 3, we solve the associated sample invariant probability vector ${\mathrm{\Pi }}_{asm}$ based on Theorem 3. For Step 2, we use the Ramaswami dual to solve the associated sample matrix $R_{asm}$. It is known that the Ramaswami dual [11] enables us to compute the matrix R for a GI/M/1-type chain with a finite phase in terms of computing the matrix G for a dual M/G/1-type Markov chain. Note that the Ramaswami dual can be modified and extended to the case of M/G/1-type and GI/M/1-type chains with a continuous phase.

5.2. Illustration with Examples

5.2.1. Example 1: An $M/G/1$-Type Chain

The Markov chain in this example is modified from Example 2 in [5] by extending the tri-diagonal structure to the more general upper-Hessenberg setting.

Denote by $S_i,i\ge 0$ the arrival times of a Poisson process with parameter $\lambda$. Let $S_0=0$. Define a sequence of i.i.d. random variables $V\left(S_n\right),n⩾0$, which are distributed with

$\mathbb{P}\left\{V\left(S_n\right)=j\right\}=p_j,j\ge -1,$

where $p_j^{\prime }s$ are non-negative constants such that ${\sum }_{j=-1}^{\infty }{p}_j=1$. We define

(15)$L\left(t\right)=\left\{\begin{array}{c}L\left(S_0\right)=L\left(0\right),\mathrm{if}0=S_0⩽t<S_1,\hfill \\ max\left\{0,L\left(S_{n-1}\right)+j\right\},\mathrm{if}V\left(S_n\right)=j,S_n⩽t<S_{n+1},\hfill \end{array}\right.$

(16)$Y\left(t\right)=\left\{\begin{array}{cc}Y\left(S_k\right)+\left(t-S_k\right),\hfill & \mathrm{if}V\left(S_k\right)=\ell ,S_k⩽t<S_{k+1},\hfill \\ t-S_k,\hfill & \mathrm{if}V\left(S_k\right)=-1,S_k⩽t<S_{k+1},\hfill \end{array}\right.$

where $n\ge 1,j\ge -1,k\ge 0,\ell \ge 0$.

Let $L_n=L\left(S_n\right)$ and $Y_n=Y\left(S_{n+1}-0\right)$, then $\left(L_n,Y_n\right)$ is a M/G/1-type chain, whose phase space is ${\mathbb{R}}_{+}$. Its transition kernels are derived as

(17)$\begin{array}{c}{A}_0(x,y)=p_{-1}\left(1-e^{-\lambda y}\right),\mathrm{for}\mathrm{all}x,y,\hfill \\ A_{j+1}(x,y)=\left\{\begin{array}{cc}0,\hfill & \mathrm{if}y<x;\hfill \\ p_j\left[1-e^{-\lambda (y-x)}\right],\hfill & \mathrm{if}y⩾x,\hfill \end{array}\right.\mathrm{if}j\ge 0.\hfill \end{array}$

and finally

(18)$B_0(x,y)=A_0(x,y)+A_1(x,y),B_i(x,y)=A_{i+1}(x,y),i\ge 1.$

The marginal invariant probability measures for the $M/G/1$-type chain $\left(L_n,Y_n\right)$ has analytical expressions, as follows:

(19)$\begin{array}{c}{}_L{\mathrm{\Pi }}_0=\frac{p_{-1}-{\sum }_{n=0}^{\infty }n{p}_n}{p_{-1}},\hfill \\ {}_L{\mathrm{\Pi }}_k=\frac{{}_L{\pi }_0\left({\sum }_{n=k}^{\infty }{p}_n\right)+{\sum }_{i=1}^{k-1}{\sum }_{j=0}^{\infty }{}_L{\pi }_i{p}_{k-i+j}}{p_{-1}},\hfill \\ {}_P\mathrm{\Pi }\left(y\right)=\left\{\begin{array}{cc}0,\hfill & \mathrm{if}y<0\hfill \\ 1-e^{-\lambda p_{-1}y},\hfill & \mathrm{if}y⩾0,\hfill \end{array}\right.\hfill \end{array}$

where ${}_L{\mathrm{\Pi }}_k$ and ${}_P\mathrm{\Pi }\left(y\right)$ are, respectively, the level and phase marginal invariant probability measures.

Take $p_{-1}=\frac{1}{2}$, $p_k={\left(\frac{1}{3}\right)}^{k+1}$, $k\ge 0$. From (19), we can have the following exact value of the marginal level stationary probabilities

${}_L{\mathrm{\Pi }}_0=\frac{1}{2},{}_L{\mathrm{\Pi }}_k=\frac{1}{4}{\left(\frac{2}{3}\right)}^k,k\ge 1.$

Evenly take 256 samples on $[0,45]$ as values of x and 256 samples on $[0,50]$ as values of y, and choose the Haar wavelet for the wavelet transform. Figure 1 presents the numerical solutions of kernel functions $G(x,y)$. The marginal distributions are obtained numerically based on the $G_{asm}$ solved by functional iteration. The numerical solutions for level marginal distribution and phase marginal distribution are, respectively, shown in Figure 2 and Figure 3, together with the corresponding analytical solutions. For each method used to derive $G_{asm}$, we calculate its mean absolute error defined as $\frac{1}{K+1}{\sum }_{k=0}^K{|}_L{\mathrm{\Pi }}_k{-}_L{\widehat{\mathrm{\Pi }}}_k|$, where $K=500$ and ${}_L{\widehat{\mathrm{\Pi }}}_k$ is the numerical solution of the marginal level stationary probability at level k. (We take $K=500$ because values of ${}_L{\mathrm{\Pi }}_k$ when $k>500$ are small enough to be considered as negligible.) In this example, different methods of solving matrix G lead to the same numerical solutions of level and phase marginal distributions. According to Figure 4, performances of various methods are similar in the sense of accuracy and computational time.

5.2.2. Example 2: A $GI/M/1$-Type Chain

Consider a first-come first-served single server $GI/G/1$ queuing system, which was considered by [1] for the theoretical analysis of the invariant probability distribution. Here, we consider the computational issue. In this $GI/G/1$ queue, the service times and interarrival times are distributed with general distribution functions $S\left(x\right)$ and $F\left(x\right)$. We assume that both the mean arrival interval $\lambda ={\int }_0^{\infty }tdF\left(t\right)$ and the mean service time ${\mu }_0:={\int }_0^{\infty }tdS\left(t\right)$ are finite.

Let $L_n$ be the number of customers right before the arrival time of the nth customer and let $X_n$ be the remaining service time just after the nth arrival. Let $Z_n$ be the departure time of the nth customer, and let $X\left(t\right)$ be the remaining service time at time t of a customer who is receiving the service. Write

$D_n^t(x,y)=P\{Z_n\le t<Z_{n+1},X\left(t\right)\le y|X\left(0\right)=x\}.$

Then, $(L_n,X_n),n\ge 1$ is a GI/M/1-type Markov chain with discrete levels and continuous phases, whose transition probabilities are given by (1) with (see [1])

(20)$A_n(x,y)={\int }_0^{\infty }{D}_n^t(x,y)dF\left(t\right),$

(21)$B_n(x,y)=\left(\sum _{j=n+1}^{\infty }{A}_j(x,\infty )\right)S\left(y\right).$

From [1], we know that if $\lambda >{\mu }_0$, then $(L_n,X_n)$ has an invariant probability measure ${\mathrm{\Pi }}_j$ with

${\mathrm{\Pi }}_k(·)=d{\int }_0^{\infty }dS\left(x\right)R^{\left(k\right)}(x,·)$

where the constant d is given by

$d=1+\left({\int }_0^{\infty }\left[\sum _{n=0}^{\infty }{F}^{n*}\left(x\right)\right]dS\left(x\right)\right)\left[exp\sum _{n=1}^{\infty }\left(1-{\int }_0^{\infty }[1-F^{n*}\left(x\right)]d{S}^{*n}\left(x\right)\right)/n\right].$

To illustrate our algorithm, we would like to compare numerical solutions for the level of marginal distribution with its analytical value. For numerical calculation, let $F\left(t\right)$ be uniformly distributed in the interval $(0,1]$, i.e., $F\left(t\right)\sim U(0,1)$, and let the service time be exponentially distributed with parameter $\mu$, i.e., $S\left(t\right)\sim exp\left(\mu \right)$. Then we have, for $n\ge 1$

$\begin{array}{c}{D}_n^t(x,y)=\left\{\begin{array}{cc}\frac{{\mu }^{n-1}{(t-x)}^{n-1}}{(n-1)!}(e^{-\mu (t-x)}-e^{-\mu (t-x+y)}),\hfill & t\ge x,\\ 0,\hfill & t<x,\end{array}\right.\end{array}$

and, for $n=0$

$D_0^t(x,y)=I_{[0,y]}(x-t).$

The kernels $A_n(x,y)$ and $B_n(x,y)$ are calculated as follows. For $n\ge 1$, by (20)

$\begin{array}{ccc}\hfill A_n(x,y)& =& {\int }_x^1\frac{{\mu }^{n-1}{(t-x)}^{n-1}}{(n-1)!}\left(e^{-\mu (t-x)}-e^{-\mu (t-x+y)}\right)dt\hfill \\ & =& \frac{1-e^{-\mu y}}{\mu }{\int }_0^{1-x}\frac{{\mu }^n}{(n-1)!}{t}^{n-1}{e}^{-\mu t}dt.\hfill \end{array}$

For $0\le x\le 1,y\ge x$

$A_0(x,y)={\int }_{x-y}^{x}1dt={\int }_0^{x}1dt=x,$

and for $0\le x\le 1,y<x$

$A_0(x,y)={\int }_{x-y}^{x}1dt=y.$

For $n\ge 0$, by (21) we have

$\begin{array}{c}\hfill B_n(x,y)=(1-e^{-\mu y})·{\int }_0^{1-x}\left[1-\sum _{j=0}^{n-1}{e}^{-\mu t}\frac{{\left(\mu t\right)}^j}{j!}\right]dt.\end{array}$

Since $\lambda =1/2$ and $E\left[S\right]={\mu }_0=\frac{1}{\mu }$, the queuing system is stable if

$\rho =\frac{1/\lambda }{1/E\left[S\right]}=\frac{E\left[S\right]}{\lambda }=\frac{2}{\mu }<1.$

It is well known that the level marginal invariant probability distribution is

${}_L{\mathrm{\Pi }}_j=c^j(1-c),$

where c is the solution of ${\int }_0^{\infty }{e}^{-\mu t(1-c)}dF\left(t\right)=c$ on the interval $[0,1]$. Since $F\sim U(0,1)$, then

$e^{(c-1)\mu }-1=c(c-1)\mu .$

For numerical experiments, we take $\mu =4.7$. The constant c is solved to be approximately 0.2885. The exact marginal level distribution ${}_L\mathrm{\Pi }$ can be obtained. On the other hand, we can perform the numerical algorithm in the previous section to approximate ${}_L\mathrm{\Pi }$. We may then compare the numerical results to the analytical results, and provide a verification of the algorithm afterward. Here, we do not consider the marginal phase distribution, since its closed form cannot be obtained.

Evenly take 256 samples on $[0,1]$ as values of x and 256 samples on $[0,1.5]$ as values of y, and choose the Haar wavelet for a wavelet transform. The numerical solutions of kernel functions $G(x,y)$ and $R(x,y)$ are shown in Figure 5. The level marginal distribution of this queuing system could be computed by a previous algorithm. We show the numerical solutions using $G_{asm}$ solved by functional iteration, together with the analytical solutions in Figure 6. Among the five numerical methods mentioned in (ii) of Remark 1, we note that the method of the invariant subspace does not work during the run of the algorithm, which may be caused by the fact that some matrices are not invertible. With numerical solutions of marginal invariant probability distributions for the level and phase, we can further estimate the mean and variance of the queue length and the remaining service time, which are listed in Table 1. This implies practical uses of invariant probability distributions.

Since we are not able to obtain analytical solutions of marginal invariant probability distributions for the phase, only the analytical mean and variance for the level are presented in Table 1. When using the Ramaswami reduction, the mean queue length is the most accurate among four methods, but the variance of the queue length is not close to the analytical variance. However, the mean and variance solved using the functional iteration are both relatively accurate.

We compare properties of the other four methods according to their mean absolute errors and speeds of computation. From Figure 7 and Table 1, the functional iteration performs the best among all the four methods, since it is the fastest and also the most accurate. When we raise the sample size from 256 to 512, it takes 20.98 s to solve $G_{asm}$ using a functional iteration. The computational times of the cyclic reduction, Newton iteration and Ramaswami reduction are 116.35 s, 1997.53 s and 2191.45 s, respectively. The differences between the mean absolute errors of numerical solutions and the errors when using a sample size of 256, however, are only around ${10}^{-4}$.

6. Conclusions

For invariant probability measures of Harris positive recurrent $GI/M/1$-type or $M/G/1$-type Markov chains with discrete levels and a general phase set, we establish wavelet-based computational frameworks in this paper. A theoretical analysis framework is also established for $M/G/1$-type Markov chains. These results extend the known findings in [4,5] for QBD processes to the current more general block-structured Markov chains. Numerical experiments support the effectiveness of our numerical algorithms based on DTWC. An interesting observation in Example 2 is that among the adopted five algorithms for G-matrix, the functional iteration performs the best, but the invariant subspace may fail.

For future research, it is interesting to consider block-structured continuous-time Markov processes with discrete levels and continuous phases. In this case, the processes should be presented in terms of the extended generators. It is expected that the research is more challenging when setting up these models and preforming the theoretical and numerical analysis of their invariant probability measures.

Author Contributions

Methodology, Y.L.; Software, S.J.; Writing—original draft, S.J., Y.L. and N.L.; Writing—review and editing, Y.L and N.L.; Visualization, S.J. and N.L.; Funding acquisition, Y.L. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

No data is used in this study.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. Numerical solution of kernel function [Forumla omitted. See PDF.].

Figure 2. Level marginal invariant probability distribution [Forumla omitted. See PDF.].

Figure 3. Phase marginal invariant probability distribution [Forumla omitted. See PDF.].

Figure 4. Difference among methods on level marginal invariant probability distribution [Forumla omitted. See PDF.]. The legend of the first plot includes mean absolute errors, and the legend of the second plot includes computational times. (FI: functional iteration, CR: cyclic reduction, NI: Newton iteration, RR: Ramaswami Reduction and IS: invariant subspace).

Figure 5. Numerical solutions of kernel functions. The right picture is about the kernel [Forumla omitted. See PDF.], and the left is about its dual kernel [Forumla omitted. See PDF.].

Figure 6. Level marginal invariant probability distribution [Forumla omitted. See PDF.].

Figure 7. Difference among methods on level marginal invariant probability distribution [Forumla omitted. See PDF.]. The legend of the first plot includes mean absolute errors, and the legend of the second plot includes computational times. (FI: functional iteration, CR: cyclic reduction, NI: Newton iteration and RR: Ramaswami Reduction).

Table 1

Mean and variance of level and phase.

	Queue Length (Level)		Remaining Service Time (Phase)
Method	Mean	Variance	Mean	Variance
FI	0.3415	0.3680	0.3125	0.1127
CR	0.6964	1.3918	0.3345	0.1255
NI	0.6964	1.3918	0.3345	0.1255
RR	0.3911	0.9596	0.4385	0.1641
Analytical	0.4054	0.5698	-	-

Appendix A

Table A1

Summary of frequently used notations.

Notation	Description
E	State space of a Markov chain
${\ell }_i$	The i level set of Markov chain $\{(L_n,X_n):n\in \mathbb{N}\}$
${\tau }_{{\ell }_i}$	The first return time to the level set ${\ell }_i$
$P_{GI}(·,·)$	Transition kernel matrix of $GI/M/1$-type Markov chain
$P_M(·,·)$	Transition kernel matrix of $M/G/1$-type Markov chain
$P^A(·,·)$	One-step transition kernel of a censored Markov chain on set A
$\mathrm{\Pi }$	Invariant probability measure
$\overline{U}$	Associated matrix of kernel function $U(x,y)$
${\sum }_H$	Set of kernel functions having a density function equaling to zeros outside of $H\times H$

References

1. Tweedie, R.L. Operator-geometric stationary distributions for Markov chains, with application to queueing models. Adv. Appl. Probab.; 1982; 14, pp. 368-391. [DOI: https://dx.doi.org/10.2307/1426527]

2. Neuts, M.F. Matrix-Geometric Solutions in Stochastic Models: An Algorithmic Approach; Johns Hopkins University Press: Baltimore, MD, USA, 1981.

3. Breuer, L. Transient and stationary distributions for the GI/G/k Queue with Lebesgue-dominated inter-arrival time distribution. Queueing Syst.; 2003; 45, pp. 47-57. [DOI: https://dx.doi.org/10.1023/A:1025643801208]

4. Nielsen, B.F.; Ramaswami, V. A Computational Framework for a Quasi Birth and Death Process with a Continuous Phase Variable; Ramaswami, V.; Wirth, P. ITC 15 Elsevier: Amsterdam, The Netherlands, 1997; pp. 477-486.

5. Jiang, S.; Latouche, G.; Liu, Y. Wavelet transform for quasi-birth-death process with a continuous phase set. Appl. Math. Comput.; 2015; 252, pp. 354-376. [DOI: https://dx.doi.org/10.1016/j.amc.2014.12.023]

6. Neuts, M.F. Structured Stochastic Matrices of M/G/1 type and Their Applications; Marcel Dekker: New York, NY, USA, 1989.

7. Bini, D.A.; Meini, B.; Steffé, S.; Van Houdt, B. Structured Markov chains solver: Software tools. Proceedings of the SMCTOOLS; Pisa, Italy, 10 October 2006.

8. Bini, D.A.; Meini, B.; Steffé, S.; Van Houdt, B. Structured Markov chains solver: Algorithms. Proceedings of the SMCTOOLS; Pisa, Italy, 10 October 2006.

9. Bright, L.W. Matrix-Analytic Methods in Applied Probability; The University of Adelaide: Adelaide, Australia, 1996.

10. Ramaswami, V. Nonlinear matrix equations in applied probability-solution techniques and open problems. Siam Rev.; 1988; 30, pp. 256-263. [DOI: https://dx.doi.org/10.1137/1030046]

11. Ramaswami, V. A duality theorem for the matrix paradigms in aueueing theory. Commun. Stat. Stoch. Model.; 1990; 6, pp. 151-161. [DOI: https://dx.doi.org/10.1080/15326349908807141]

12. Taylor, P.G.; Van Houdt, B. On The dual relationship between Markov chains of GI/M/1 and M/G/1 type. Adv. Appl. Probab.; 2010; 42, pp. 210-225. [DOI: https://dx.doi.org/10.1239/aap/1269611150]

13. Bini, D.A.; Latouche, G.; Meini, B. Numerical Methods for Structured Markov Chains; Oxford University Press: New York, NY, USA, 2005.

14. Meyn, S.P.; Tweedie, R.L. Markov Chains And Stochastic Stability; 2nd ed. Cambridge University Press: Cambridge, UK, 2009.

15. Stéphane, M. A Wavelet Tour of Signal Processing; Academic Press: New York, NY, USA, 2009.

16. Daubechies, I. Orthonormal base of compactly supported wavelets. Commun. Pure Appl. Math.; 1988; 41, pp. 909-996. [DOI: https://dx.doi.org/10.1002/cpa.3160410705]

17. Favati, P.; Meini, B. Relaxed functional iteration techniques for the numerical solution of M/G/1 type Markov chains. Bit Numer. Math.; 1998; 38, pp. 510-526. [DOI: https://dx.doi.org/10.1007/BF02510257]

18. Favati, P.; Meini, B. On functional iteration methods for solving nonlinear matrix equations arising in queueing problems. IMA J. Numer. Anal.; 1999; 19, pp. 39-49. [DOI: https://dx.doi.org/10.1093/imanum/19.1.39]

19. Latouche, G. Newton’s iteration for non-linear equations in Markov chains. IMA J. Numer. Anal.; 1994; 14, pp. 583-598. [DOI: https://dx.doi.org/10.1093/imanum/14.4.583]

20. Akar, N.; Sohraby, K. An invariant subspace approach in M/G/1 and G/M/1 type Markov chains. Stoch. Model.; 1997; 13, pp. 381-416. [DOI: https://dx.doi.org/10.1080/15326349708807433]

21. Bini, D.A.; Meini, B. On Cyclic Reduction Applied to a Class of Toeplitz-like Matrices Arising in Queueing Problems, Computations with Markov Chains; Springer: Boston, MA, USA, 1995.

22. Ramaswami, V. The generality of Quasi Birth-and-Death processes. Advances in Matrix Analytic Methods for Stochastic Models; Alfa, A.S.; Chakravarthy, S.R. Notable Publications: Branchburg, NJ, USA, 1998; pp. 93-113.