Dynamic Analysis of a Unified Multivariate Counting Process and Its Asymptotic Behavior

Abstract

The class of counting processes constitutes a significant part of applied probability. The classic counting processes include Poisson processes, nonhomogeneous Poisson processes, and renewal processes. More sophisticated counting processes, including Markov renewal processes, Markov modulated Poisson processes, age-dependent counting processes, and the like, have been developed for accommodating a wider range of applications. These counting processes seem to be quite different on the surface, forcing one to understand each of them separately. The purpose of this paper is to develop a unified multivariate counting process, enabling one to express all of the above examples using its components, and to introduce new counting processes. The dynamic behavior of the unified multivariate counting process is analyzed, and its asymptotic behavior as 𝑡→∞ is established. As an application, a manufacturing system with certain maintenance policies is considered, where the optimal maintenance policy for minimizing the total cost is obtained numerically.

1. Introduction

A stochastic process {𝑁(𝑡)∶𝑡≥0} is called a counting process when 𝑁(𝑡) is nonnegative, right continuous and monotone nondecreasing with 𝑁(0)=0. The classic counting processes of importance include a Poisson process, a nonhomogeneous Poisson process (NHPP) and a renewal process. More sophisticated counting processes have been developed in order to accommodate a wider range of applications. A Markov renewal process, for example, extends an ordinary renewal process in that the interarrival time between two successive arrivals has a probability distribution depending on the state transition of the underlying Markov chain, see, for example, Pyke [1, 2]. In Masuda and Sumita [3], the number of entries of a semi-Markov process into a subset of the state space is analyzed, while a Markov-modulated Poisson process (MMPP) is introduced by Neuts [4] where jumps of the MMPP occur according to a Poisson process with intensity 𝜆_𝑖 whenever the underlying Markov-chain is in state 𝑖.

The MMPP is generalized subsequently into several different directions. Lucantoni et al. [5] develop a Markovian arrival process (MAP), where a Markov-chain defined on 𝒥=𝐺∪𝐵 with 𝐺∩𝐵=𝜙,𝐺≠𝜙, and 𝐵≠𝜙 is replaced to a state 𝑖∈𝐺 as soon as it enters a state 𝑗∈𝐵 with probability 𝑝̃_𝑗𝑖 and the counting process describes the number of such replacements occurred in [0,𝑡]. If we set 𝐺={0_𝐺,1_𝐺,…,𝐽_𝐺} and 𝐵={0_𝐵,1_𝐵,…,𝐽_𝐵}, and an absorbing state 𝑗_𝐵∈𝐵 can be reached only from its counter part 𝑗_𝐺∈𝐺 with instantaneous replacement back to 𝑗_𝐺 itself, the resulting MAP becomes an MMPP. Another direction for generalizing the MMPP is to replace the underlying Markov-chain by a semi-Markov process, which we call a semi-Markov-modulated Poisson process (SMMPP). To the best knowledge of the authors, the SMMPP is first addressed by Dshalalow [6], where a systematic approach for dealing with modulated random measures is proposed in a more abstract way. Such modulated random measures include the MMPP and the SMMPP as special cases. An application of the SMMPP to queueing theory is discussed in Dshalalow and Russell [7]. More recently, in a series of papers by Agarwal et al. [8] and Dshalalow [9, 10], the original approach above has been further extended and refined. The SMMPP is studied in detail independently by Özekici and Soyer [11].

In the SMMPP, the counting process under consideration is modulated according to state transitions of the underlying semi-Markov process. A further generalization may be possible by considering a counting process whose arrival intensity depends on not only the current state of the semi-Markov process but also the current age of the process. This line of research is originated by Sumita and Shanthikumar [12] where an age-dependent counting process generated from a renewal process is studied. Here, items arrive according to an NHPP which is interrupted and reset at random epochs governed by a renewal process.

All of these counting processes discussed above seem to be quite different on the surface, forcing one to understand each of them separately. The purpose of this paper is to develop a unified multivariate counting process which would contain all of the above examples as special cases. In this regard, we consider a system where items arrive according to an NHPP. This arrival stream is interrupted from time to time where the interruptions are governed by a finite semi-Markov process 𝐽(𝑡) on 𝒥={0,1,2,…,𝐽}. Whenever a state transition of the semi-Markov process occurs from 𝑖 to 𝑗, the intensity function of the NHPP is switched from 𝜆_𝑖(𝑥) to 𝜆_𝑗(𝑥) with an initial value reset to 𝜆_𝑗(0). In other words, the arrivals of items are generated by the NHPP with 𝜆_𝑖(𝑥) when the semi-Markov process is in state 𝑖 with 𝑥 denoting the time since the last entry into state 𝑖. Of particular interest in analysis of such systems are the multivariate counting processes 𝑀(𝑡)=[𝑀_𝑖(𝑡)]_𝑖∈𝒥 and 𝑁(𝑡)=[𝑁_𝑖𝑗(𝑡)]_{𝑖,𝑗∈𝒥} where 𝑀_𝑖(𝑡) counts the cumulative number of items that have arrived in [0,𝑡] while the semi-Markov process is in state 𝑖 and 𝑁_𝑖𝑗(𝑡) represents the cumulative number of the state transitions of the semi-Markov process from 𝑖 to 𝑗 in [0,𝑡]. The joint multivariate counting process [𝑀(𝑡),𝑁(𝑡)] enables one to unify many existing counting processes in that they can be derived in terms of the components of [𝑀(𝑡),𝑁(𝑡)]. Because of this reason, hereafter, we call [𝑀(𝑡),𝑁(𝑡)] the unified multivariate counting process. The dynamic behavior of [𝑀(𝑡),𝑁(𝑡)] is captured through analysis of the underlying Laplace transform generating functions, yielding asymptotic expansions of the means and the variances of the partial sums of its components as 𝑡→∞.

Applications of the unified multivariate counting process can be found, for example, in modern communication networks. One may consider a high-speed communication link for transmitting video signals between two locations. Video sequences are transmitted as streams of binary data that vary over time in traffic intensity according to the level of movement, the frequency of scene changes, and the level of transmission quality. Consequently, efficient transmission of video traffic can be achieved through variable bit rate coding. In this coding scheme, data packets are not generated at a constant rate from the original sequence, but rather at varying rates. By doing so, one achieves less fluctuation in transmission quality level and, at the same time, transmission capacity can be freed up whenever possible. As in Maglaris et al. [13], such a mechanism may be implemented by using multimode encoders where each mode reflects a certain level of data compression, and the change between modes is governed by the underlying video sequence according to buffer occupancy levels. A system of this sort can be described in the above framework with 𝑀(𝑡)=∑_𝑖∈𝒥𝑀_𝑖(𝑡) representing the number of packet arrivals at the origination site and 𝑁_𝑖𝑗(𝑡) describing the number of the encoder changes in [0,𝑡]. The state of the underlying semi-Markov process at time 𝑡 then corresponds to the current mode of the encoder. Other types of applications include system reliability models where the semi-Markov process describes the status of the system under consideration while the interruptions correspond to system failures and replacements. A cost function associated with such a system may then be constructed from the unified multivariate counting processes [𝑀(𝑡),𝑁(𝑡)]. In this paper, a manufacturing system with certain maintenance policies is considered, where the unified multivariate counting process enables one to determine numerically the optimal maintenance policy for minimizing the total cost.

The structure of this paper is as follows. In Section 2, key transform results of various existing counting processes are summarized. Detailed description of the unified multivariate counting process is provided in Section 3 and its dynamic behavior is analyzed in Section 4 by examining the probabilistic flow of the underlying stochastic processes and deriving transform results involving Laplace transform generating functions. Section 5 is devoted to derivation of the existing counting processes of Section 2 in terms of the components of the unified multivariate counting process. Asymptotic analysis is provided in Section 6, yielding asymptotic expansions of the means and the variances of the partial sums of its components. An application is discussed in Section 7, where a manufacturing system with certain maintenance policies is considered and the optimal maintenance policy for minimizing the total cost is obtained numerically. Some concluding remarks are given in Section 8.

Throughout the paper, matrices and vectors are indicated by double underlines (𝐴,𝑏, etc.) and underlines (𝑋,𝑦, etc.) respectively. The vector with all components equal to 1 is denoted by 1 and the identity matrix is written as 𝐼. A submatrix of 𝑎 is defined by 𝑎_𝐺𝐵=[𝑎_𝑖𝑗]_{𝑖∈𝐺,𝑗∈𝐵}.

2. Various Counting Processes of Interest

In this section, we summarize key transform results of various counting processes of interest, which can be expressed in terms of the components of the unified multivariate counting process proposed in this paper as we will see. We begin the discussion with one of the most classical arrival processes, the Poisson process.

2.1. Poisson Process

Poisson process of intensity 𝜆 is characterized by a sequence of independently and identically distributed (i.i.d.) exponential random variables (𝑋_𝑗)^∞𝑗=1 with common probability density function (p.d.f.) 𝑓_𝑋(𝑥)=𝜆𝑒^−𝜆𝑥. Let 𝑆_𝑛=∑^𝑛𝑗=1𝑋_𝑗. Then, the associated Poisson process {𝑁(𝑡)∶𝑡≥0} is defined as a counting process satisfying 𝑁(𝑡)=𝑛⟺𝑆_𝑛≤𝑡<𝑆_𝑛+1.(2.1) If a system has an exponential lifetime of mean 𝜆⁻¹ and is renewed instantaneously upon failure, 𝑋_𝑗 represents the lifetime of the 𝑗th renewal cycle. The Poisson process {𝑁(𝑡)∶𝑡≥0} then counts the number of failures that have occurred by time 𝑡.

Let 𝑝_𝑛(𝑡)=P[𝑁(𝑡)=𝑛∣𝑁(0)=0] and define the probability generating function (p.g.f.) 𝜋(𝑣,𝑡) by 𝜋(𝑣,𝑡)=E𝑣^𝑁(𝑡)=^∞_𝑛=0𝑝_𝑛(𝑡)𝑣^𝑛.(2.2) It can be seen, see, for example, Karlin and Taylor [14], that 𝑑𝑑𝑡𝑝^𝑛(𝑡)=−𝜆𝑝_𝑛(𝑡)+𝜆𝑝_𝑛−1(𝑡)(2.3) where 𝑝_𝑛(𝑡)=0 for 𝑛<0. Multiplying 𝑣^𝑛 on both sides of (2.3) and summing from 0 to ∞, one then finds that 𝜕𝜕𝑡𝜋(𝑣,𝑡)=−𝜆(1−𝑣)𝜋(𝑣,𝑡).(2.4) Since 𝑝_𝑛(0)=𝛿_{𝑛=0} where 𝛿_{𝑃}=1 if statement 𝑃 is true and 𝛿_{𝑃}=0 otherwise, one has 𝜋(𝑣,0)=1. Equation (2.4) can then be solved as 𝜋(𝑣,𝑡)=𝑒^{−𝜆𝑡(1−𝑣)};𝑝_𝑛(𝑡)=𝑒^−𝜆𝑡(𝜆𝑡)^𝑛𝑛!.(2.5)

2.2. Nonhomogeneous Poisson Process (NHPP)

An NHPP {𝑀(𝑡)∶𝑡≥0} differs from a Poisson process in that the failure intensity of the system is given as a function of time 𝑡. Accordingly, (2.3) should be rewritten as 𝑑𝑑𝑡𝑝^𝑚(𝑡)=−𝜆(𝑡)𝑝_𝑚(𝑡)+𝜆(𝑡)𝑝_𝑚−1(𝑡).(2.6) By taking the generating function of (2.6), one finds that 𝜕𝜕𝑡𝜋(𝑢,𝑡)=−𝜆(𝑡)(1−𝑢)𝜋(𝑢,𝑡).(2.7) With 𝐿(𝑡)=∫^𝑡0𝜆(𝑦)𝑑𝑦, this equation can be solved as 𝜋(𝑢,𝑡)=𝑒^{−𝐿(𝑡)(1−𝑢)};𝑝_𝑚(𝑡)=𝑒^{−𝐿(𝑡)}𝐿(𝑡)^𝑚𝑚!.(2.8) The reader is referred to Ross [15] for further discussions of NHPPs.

2.3. Markov-Modulated Poisson Process (MMPP)

Let {𝐽(𝑡)∶𝑡≥0} be a Markov-chain in continuous time on 𝒥={0,…,𝐽} governed by a transition rate matrix 𝜈=[𝜈_𝑖𝑗]. Let 𝜆^⊤=[𝜆₀,…,𝜆_𝐽] and define the associated diagonal matrix 𝜆_𝐷=[𝛿_{𝑖=𝑗}𝜆_𝑖]. An MMPP {𝑀(𝑡)∶𝑡≥0} characterized by (𝜈,𝜆_𝐷) is a pure jump process where jumps of 𝑀(𝑡) occur according to a Poisson process with intensity 𝜆_𝑖 whenever the Markov-chain 𝐽(𝑡) is in state 𝑖.

Let 𝜈_𝑖=∑_𝑗∈𝒥𝜈_𝑖𝑗 and define 𝜈_𝐷=[𝛿_{𝑖=𝑗}𝜈_𝑖]. The infinitesimal generator 𝑄 associated with the Markov-chain 𝐽(𝑡) is then given by 𝑄=−𝜈_𝐷+𝜈.(2.9) For 𝑖,𝑗∈𝒥, let 𝑝𝑝(𝑘,𝑡)=𝑝_𝑖𝑗(𝑘,𝑡)_𝑖𝑗;(𝑘,𝑡)=P[𝑀(𝑡)=k,𝐽(𝑡)=𝑗∣𝐽(0)=𝑖,𝑀(0)=0],(2.10) and define the associated matrix generating function 𝜋(𝑢,𝑡) by 𝜋(𝑢,𝑡)=^∞_𝑘=00𝑥0200𝑑𝑝(𝑘,𝑡)𝑢^𝑘.(2.11) It can be seen that 𝜕𝜕𝑡𝑝^𝑖𝑗(𝑘,𝑡)=−𝜆_𝑗+𝜈_𝑗𝑝_𝑖𝑗(𝑘,𝑡)+_𝑟∈𝒥0𝑥0200𝑑𝑝_𝑖𝑟(𝑘,𝑡)𝜈_𝑟𝑗+𝜆_𝑗𝑝_𝑖𝑗(𝑘−1,𝑡).(2.12) In matrix notation, this can be rewritten as 𝜕𝜕𝑡𝑝(𝑘,𝑡)=−𝑝(𝑘,𝑡)𝜆_𝐷+𝜈_𝐷−𝜈+𝑝(𝑘−1,𝑡)𝜆_𝐷.(2.13) By taking the generating function of (2.13) together with (2.9), one sees that 𝜕𝜕𝑡𝜋(𝑢,𝑡)=𝜋(𝑢,𝑡)𝑄−(1−𝑢)𝜆_𝐷.(2.14) Since 𝑀(0)=0, one has 𝜋(𝑢,0)=𝐼, where 𝐼=[𝛿_{𝑖=𝑗}] is the identity matrix, so that the above differential equation can be solved as 𝜋(𝑢,𝑡)=𝑒^{𝑄−(1−𝑢)𝜆𝐷𝑡}=^∞_𝑘=0𝑡^𝑘𝑘!𝑄−(1−𝑢)𝜆_𝐷^𝑘,(2.15) where 𝐴=^0def𝐼 for any square matrix 𝐴. It should be noted that 𝜋(1,𝑡)=𝑒^𝑄𝑡, which is the transition probability matrix of 𝐽(𝑡) as it should be. By taking the Laplace transform of both sides of (2.15), 𝜋(𝑢,𝑠)=∫^∞0𝑒^−𝑠𝑡𝜋(𝑢,𝑡) is given by 𝜋(𝑢,𝑠)=𝑠𝐼−𝑄+(1−𝑢)𝜆_𝐷⁻¹.(2.16)

In general, the interarrival times generated by an MMPP are not independent nor identically distributed. In multimedia computer and communication networks, data packets are mingled together with voice and image packets generated from analogue sources. Since arrival patterns of such packets differ from each other, MMPPs have provided useful means to model arrival processes of packets in multimedia computer and communication networks, see, for example, Heffes and Lucantoni [16] and Sriram and Whitt [17]. Burman and Smith [18], and Knessl et al. [19] studied a single server queuing system with an MMPP arrival process and general i.i.d. service times. Neuts et al. [20] established characterization theorems for an MMPP to be a renewal process in terms of lumpability of the underlying Markov-chain 𝐽(𝑡). The reader is referred to Neuts [4] for further discussions of MMPP.

An MMPP can be extended by replacing the underlying Markov-chain in continuous time by a semi-Markov process as discussed in Section 1. This process is denoted by SMMPP. To the best knowledge of the authors, the SMMPP is first addressed by Dshalalow [6], where a systematic approach for dealing with modulated random measures is proposed in a more abstract way. Such modulated random measures include the MMPP and the SMMPP as special cases. An application of the SMMPP to queueing theory is addressed in Dshalalow and Russell [7]. More recently, the original approach above has been further extended and refined in a series of papers by Agarwal et al. [8] and Dshalalow [9, 10]. The SMMPP is studied in detail independently by Özekici and Soyer [11] including transient characterizations and ergodic analysis. Both MMPP and SMMPP will be proven to be expressible in terms of the components of the unified multivariate counting process proposed in this paper.

2.4. Renewal Process

Renewal processes can be considered as a generalization of Poisson processes in that a sequence of i.i.d. exponential random variables are replaced by that of any i.i.d. nonnegative random variables with common distribution function 𝐴(𝑥). The resulting counting process {𝑁(𝑡)∶𝑡≥0} is still characterized by (2.1). Let 𝑝_𝑛(𝑡)=P[𝑁(𝑡)=𝑛∣𝑁(0)=0] as before. One then sees that 𝑝_𝑛(𝑡)=𝐴^(𝑛)(𝑡)−𝐴^(𝑛+1)(𝑡),(2.17) where 𝐴^(𝑛)(𝑡) denotes the 𝑛-fold convolution of 𝐴(𝑥) with itself, that is, 𝐴^(𝑛+1)(𝑡)=∫^𝑡0𝐴^(𝑛)(𝑡−𝑥)𝑑𝐴(𝑥) and 𝐴⁽⁰⁾(𝑡)=𝑈(𝑡) which is the step function defined as 𝑈(𝑡)=1 for 𝑡≥0 and 𝑈(𝑡)=0 else.

Let 𝜋_𝑛(𝑠)=∫^∞0𝑒^−𝑠𝑡𝑝_𝑛(𝑡)𝑑𝑡 and 𝛼(𝑠)=∫^∞0𝑒^−𝑠𝑡𝑑𝐴(𝑡). By taking the Laplace transform of both sides of (2.17), it follows that 𝜋_𝑛(𝑠)=1−𝛼(𝑠)𝑠𝛼(𝑠)^𝑛.(2.18) By taking the generating function of the above equation with 𝜋(𝑣,𝑠)=∑^∞𝑛=0𝜋_𝑛(𝑠)𝑣^𝑛, one has 𝜋(𝑣,𝑠)=1−𝛼(𝑠)𝑠1−𝑣1𝛼(𝑠).(2.19) The reader is referred to Cox [21], or Karlin and Taylor [14] for further discussions of renewal processes.

2.5. Markov Renewal Process (MRP)

An MRP is an extension of an ordinary renewal process in that, in the interval [0,𝑡), the former describes the recurrence statistics for intermingling classes of epochs of an underlying semi-Markov process, whereas the latter counts the number of recurrences for a single recurrent class of epochs. More specifically, let {𝐽(𝑡)∶𝑡≥0} be a semi-Markov process on 𝒥={0,…,𝐽} governed by a matrix p.d.f. 𝑎(𝑥) where 𝑎(𝑥)≥0 and ∫^∞0𝑎(𝑥)𝑑𝑥=𝐴₀ is a stochastic matrix which is assumed to be ergodic. Let 𝜀(ℓ) be a recurrent class consisting of the entries of the semi-Markov process to state ℓ for ℓ∈𝒥, and define 𝑁_ℓ𝑟(𝑡) to be a counting process describing the number of recurrences for 𝜀(𝑟) given that there was an epoch of 𝜀(ℓ) at time 𝑡=0. Then 𝑁_ℓ(𝑡)=[𝑁_ℓ0(𝑡),…,𝑁_ℓ𝐽(𝑡)] is called an MRP.

The study of MRPs can be traced back to early 1960s represented by the two original papers by Pyke [1, 2], followed by Keilson [22, 23], Keilson and Wishart [24, 25], Çinlar [26, 27] and McLean and Neuts [28]. Since then, the area attracted many researchers and a survey paper by Çinlar [29] in 1975 already included more than 70 leading references. The study has been largely focused on the matrix renewal function 𝐻(𝑡)=[𝐻_ℓ𝑟(𝑡)] with 𝐻_ℓ𝑟(𝑡)=E[𝑁_ℓ𝑟(𝑡)], the associated matrix renewal density, and the limit theorems. For example, one has the following result concerning the Laplace transform of 𝐻(𝑡) by Keilson [23]: ℒ𝐻(𝑡)=1𝑠𝛼(𝑠)𝐼−𝛼(𝑠)⁻¹,(2.20) where 𝛼(𝑠) is the Laplace transform of 𝑎(𝑡). The unified multivariate counting process of this paper contains an MRP as a special case and provides more information based on dynamic analysis of the underlying probabilistic flows.

2.6. Number of Entries of a Semi-Markov Process into a Subset of the State Space (NESMPS)

Another type of counting processes associated with a semi-Markov process on 𝒥={0,…,𝐽} governed by a matrix p.d.f. 𝑎(𝑥) is studied in Masuda and Sumita [3], where the state space 𝒥 is decomposed into a set of good states 𝐺(≠𝜙) and a set of bad states 𝐵(≠𝜙) satisfying 𝒥=𝐺∪𝐵 and 𝐺∩𝐵=𝜙. The counting process 𝑁_𝐺𝐵(𝑡) is then defined to describe the number of entries of 𝐽(𝑡) into 𝐵 by time 𝑡.

While 𝑁_𝐺𝐵(𝑡) is a special case of MRPs, the detailed analysis is provided in [3], yielding much more information. More specifically, let 𝑋(𝑡) be the age process associated with 𝐽(𝑡), that is, 𝑋(𝑡)=𝑡−sup𝜏∶𝐽(𝑡)∣^𝜏+𝜏−≠0,0<𝜏≤𝑡,(2.21) where 𝑓(𝑥)∣^𝑥+𝑥−=𝑓(𝑥+)−𝑓(𝑥−), and define 𝐹_𝑛(𝑥,𝑡)=𝐹_{𝑛∶𝑖𝑗}(𝑥,𝑡),(2.22) where 𝐹_{𝑛∶𝑖𝑗}(𝑥,𝑡)=P𝑋(𝑡)≤𝑥,𝑁_𝐺𝐵(𝑡)=𝑛,𝐽(𝑡)=𝑗∣𝑋(0)=𝑁_𝐺𝐵(0)=0,𝐽(0)=𝑖.(2.23) One then has 𝑓_𝑛(𝑥,𝑡)=𝜕𝜕𝑥𝐹^𝑛(𝑥,𝑡).(2.24) The associated matrix Laplace transform generating function can then be defined as 𝜑(𝑣,𝑤,𝑠)=^∞_𝑛=0𝑣^𝑛^∞0𝑒^{−𝑤𝑡−𝑠𝑡}𝑓_𝑛(𝑥,𝑡)𝑑𝑥𝑑𝑡.(2.25) It has been shown in Masuda and Sumita [3] that 𝜑(𝑣,𝑤,𝑠)=1𝑤+𝑠⋅𝛾₀(𝑠)𝐼−𝑣𝛽(𝑠)⁻¹𝐼−𝛼_𝐷(𝑤+𝑠).(2.26) Here, with 𝛼(𝑠)=∫^∞0𝑒^−𝑠𝑡𝑎(𝑡)𝑑𝑡 and 𝛼_𝐶𝐷(𝑠)=[𝛼_𝑖𝑗(𝑠)]_{𝑖∈𝐶,𝑗∈𝐷} for 𝐶,𝐷⊂𝒥, the following notation is employed:

𝛼(𝑠)=𝛼^𝐺𝐺𝛼(𝑠)𝛼_𝐺𝐵(𝑠)_𝐵𝐺𝛼(𝑠)𝛼_𝐵𝐵(𝑠)_𝐷;𝜒(𝑠)=𝛿_{𝑖=𝑗}𝛼_𝑖(𝑠)with𝛼_𝑖(𝑠)=∑_𝑗∈𝒥𝛼_𝑖𝑗(𝑠);(2.27)_𝐺(𝑠)=𝛽𝐼_𝐺𝐺−𝛼_𝐺𝐺(𝑠)⁻¹;𝜒_𝐵(𝑠)=𝐼_𝐵𝐵−𝛼_𝐵𝐵(𝑠)⁻¹;(2.28)(𝑠)=0^𝐵𝐵0_𝐵𝐺𝛾𝛼_𝐺𝐵𝜒_𝐵(𝑠)𝛼_𝐺𝐵𝜒_𝐵(𝑠)𝛼_𝐵𝐺𝜒_𝐺(𝑠);(2.29)₀(𝑠)=𝛼^𝐵𝐵0𝜒_𝐵(𝑠)𝜒_𝐵(𝑠)𝛼_𝐵𝐺𝜒_𝐺(𝑠)_𝐺𝐵𝛼_𝐺𝐺𝜒_𝐺(𝑠)+𝐼.(2.30) As we will see, the unified multivariate counting process proposed in this paper enables one to deal with multidimensional generalization of NESMPSs as a special case.

2.7. Markovian Arrival Process (MAP)

As for Poisson processes, a renewal process requires interarrival times to form a sequence of i.i.d. nonnegative random variables. As we have seen, a class of MMPPs enables one to avoid this requirement by introducing different Poisson arrival rates depending on the state of the underlying Markov chain. An alternative way to avoid this i.i.d. requirement is to adapt a class of MAPs, originally introduced by Lucantoni et al. [5]. We discuss here a slightly generalized version of MAPs in that a set of absorbing states is not necessarily a singleton set.

Let {𝐽^∗(𝑡)∶𝑡≥0} be an absorbing Markov-chain on 𝒥^∗=𝐺∪𝐵 with 𝐺≠𝜙,𝐵≠𝜙 and 𝐺∩𝐵=𝜙, where all states in 𝐵 are absorbing. Without loss of generality, we assume that 𝐺={0,…,𝑚}, and 𝐵={𝑚+1,…,𝑚+𝐾}. For notational convenience, the following transition rate matrices are introduced.

𝜈^∗𝐺𝐺=𝜈_𝑖𝑗_{𝑖,𝑗∈𝐺};𝜈^∗𝐺𝐵=𝜈_𝑖𝑟_{𝑖∈𝐺,𝑟∈𝐵}.(2.31) The entire transition rate matrix 𝜈^∗ governing 𝐽^∗(𝑡) is then given by 𝜈^∗=𝜈^∗𝐺𝐺𝜈^∗𝐺𝐵00.(2.32)

A replacement Markov-chain {𝐽(𝑡)∶𝑡≥0} on 𝐺 is now generated from {𝐽^∗(𝑡)∶𝑡≥0}. Starting from a state in 𝐺, the process 𝐽(𝑡) coincides with 𝐽^∗(𝑡) within 𝐺. As soon as 𝐽^∗(𝑡) reaches state 𝑟∈𝐵, it is instantaneously replaced at state 𝑗∈𝐺 with probability 𝑝̃_𝑟𝑗 and the process continues. Let 𝐶=𝜈^∗𝐺𝐺,𝐷=𝜈^∗𝐺𝐵𝑝̃_𝐵𝐺,(2.33) where 𝑝̃_𝐵𝐺=[𝑝̃_𝑟𝑗]_{𝑟∈𝐵,𝑗∈𝐺}. Then the transition rate matrix 𝜈 and the infinitesimal generator 𝑄 of 𝐽(𝑡) are given as 𝜈=𝐶+𝐷;𝑄=−𝜈_𝐷+𝜈,(2.34) where 𝜈_𝐷=𝐶_𝐷+𝐷_𝐷,(2.35) with 𝐷𝐶_𝐷𝐷=𝛿_{𝑖=𝑗}𝑐_𝑖;𝑐_𝑖=_𝑗∈𝐺𝑐_𝑖𝑗,=𝛿_{𝑖=𝑗}𝑑_𝑖;𝑑_𝑖=_𝑗∈𝐺𝑑_𝑖𝑗.(2.36) Let 𝑝(𝑡) be the transition probability matrix of 𝐽(𝑡) with its Laplace transform defined by 𝜋(𝑠)=∫^∞0𝑒^−𝑠𝑡𝑝(𝑡)𝑑𝑡. From the Kolmogorov forward equation (𝑑/𝑑𝑡)𝑝(𝑡)=𝑝(𝑡)𝑄 with 𝑝(0)=𝐼, one has 𝜋(𝑠)=𝑠𝐼−𝑄⁻¹.(2.37)

Let {𝑁_MAP(𝑡)∶𝑡≥0} be the counting process keeping the record of the number of replacements in [0,𝑡) and define 𝑓_𝑖𝑗(𝑘,𝑡)=𝑃𝑁_MAP(𝑡)=𝑘,𝐽(𝑡)=𝑗∣𝐽(0)=𝑖,𝑖,𝑗∈𝐺.(2.38) By analyzing the probabilistic flow at state 𝑗 at time 𝑡+Δ, it can be seen that 𝑓_𝑖𝑗(𝑘,𝑡+Δ)=+𝑓_𝑖𝑗(_ℓ∈𝐺𝑘,𝑡)1−_ℓ∈𝐺𝑐_𝑗ℓ+𝑑_𝑗ℓΔ+_ℓ∈𝐺𝑓_𝑖ℓ(𝑘,𝑡)𝑐_ℓ𝑗Δ𝑓_𝑖ℓ(𝑘−1,𝑡)𝑑_ℓ𝑗Δ+𝑜(Δ).(2.39) It then follows that 𝜕𝜕𝑡𝑓^𝑖𝑗(𝑘,𝑡)=−𝑓_𝑖𝑗(𝑘,𝑡)_ℓ∈𝐺𝑐_𝑗ℓ+𝑑_𝑗ℓ+_ℓ∈𝐺𝑓_𝑖ℓ(𝑘,𝑡)𝑐_ℓ𝑗+_ℓ∈𝐺𝑓_𝑖ℓ(𝑘−1,𝑡)𝑑_ℓ𝑗.(2.40) In matrix notation, the above equation can be rewritten as 𝜕𝜕𝑡𝑓(𝑘,𝑡)=−𝑓(𝑘,𝑡)𝜈_𝐷+𝑓(𝑘,𝑡)𝐶+𝑓(𝑘−1,𝑡)𝐷.(2.41) We now introduce the following matrix Laplace transform generating function: 𝜑(𝑣,𝑠)=𝜑_𝑖𝑗(𝑣,𝑠);𝜑_𝑖𝑗(𝑣,𝑠)=^∞0𝑒^{∞−𝑠𝑡}_𝑘=0𝑓_𝑖𝑗(𝑘,𝑡)𝑣^𝑘𝑑𝑡.(2.42) Taking the Laplace transform with respect to 𝑡 and the generating function with respect to 𝑘 of both sides of (2.41), one has 𝜑(𝑣,𝑠)𝑠𝐼+𝜈_𝐷−𝐶−𝑣𝐷=𝐼,(2.43) which can be solved as 𝜑(𝑣,𝑠)=𝑠𝐼−𝑄+(1−𝑣)𝐷⁻¹.(2.44) We note from (2.37) and (2.44) that 𝜑(1,𝑠)=𝜋(𝑠) as it should be.

It may be worth noting that an MMPP is a special case of an MAP, which can be seen in the following manner. Let an MAP be defined on 𝐺∪𝐵 where 𝐺={0_𝐺,…,𝐽_𝐺} and 𝐵={0_𝐵,…,𝐽_𝐵}. Transitions within 𝐺 is governed by 𝜈. An entry into 𝑗_𝐵∈𝐵 is possible only from 𝑗_𝐺∈𝐺. When this occurs, the Markov process is immediately replaced at the entering state 𝑗_𝐺. The counting process for the number of such replacements then has the Laplace transform generating function 𝜑(𝑣,𝑠) given in (2.44) where 𝐷 is replaced by 𝜆_𝐷, which coincides with 𝜋(u,𝑠) of (2.16) for MMPPs, proving the claim.

2.8. Age-Dependent Counting Process Generated from a Renewal Process (ACPGRP)

An age-dependent counting process generated from a renewal process has been introduced and studied by Sumita and Shanthikumar [12], where items arrive according to an NHPP which is interrupted and reset at random epochs governed by a renewal process. More specifically, let {𝑁(𝑡)∶𝑡≥0} be a renewal process associated with a sequence of i.i.d. nonnegative random variables (𝑋_𝑗)^∞𝑗=1 with common p.d.f. 𝑎(𝑥). The age process 𝑋(𝑡) is then defined by 𝑋(𝑡)=𝑡−sup𝜏∶𝑁(𝑡)∣^𝜏+𝜏−=1,0<𝜏≤𝑡.(2.45) In other words, 𝑋(𝑡) is the elapsed time since the last renewal. We next consider an NHPP 𝑍(𝑥) governed by an intensity function 𝜆(𝑥). If we define 𝐿(𝑥)=^𝑥0𝜆(𝑦)𝑑𝑦,(2.46) one has 𝑔(𝑥,𝑘)=P[𝑍(𝑥)=𝑘]=exp(−𝐿(𝑥))𝐿(𝑥)^𝑘𝑘!,𝑘=0,1,2,….(2.47) Of interest, then, is a counting process {𝑀(𝑡)∶𝑡≥0} characterized by P[𝑀=(𝑡𝜆+(𝑥Δ))Δ−+𝑀𝑜((Δ𝑡))=,1(𝑚∣𝑀,𝑛,(𝑥𝑡))=∈𝑚𝑆,,𝑁Δ>(𝑡)0=.1,𝑋(𝑡)=𝑥](2.48) Here 𝑆=ℤ₊×ℤ₊×ℝ₊ where ℤ₊ is the set of nonnegative integers, and ℝ₊ is the set of nonnegative real numbers.

Since the counting process 𝑀(𝑡) is governed by the intensity function depending on the age process 𝑋(𝑡) of the renewal process 𝑁(𝑡), it is necessary to analyze the trivariate process [𝑀(𝑡),𝑁(𝑡),𝑋(𝑡)]. Let the multivariate transform of [𝑀(𝑡),𝑁(𝑡),𝑋(𝑡)] be defined by 𝜑(𝑢,𝑣,𝑤,𝑠)=^∞0𝑒^−𝑠𝑡E𝑢^𝑀(𝑡)𝑣^𝑁(𝑡)𝑒^{−𝑤𝑋(𝑡)}𝑑𝑡.(2.49) It has been shown in Sumita and Shanthikumar [12] that 𝜑(𝑢,𝑣,𝑤,𝑠)=𝛽^∗(𝑢,𝑤+𝑠)1−𝑣𝛽(𝑢,𝑠),(2.50) where we define, for 𝑚≥0 with 𝐴(𝑡)=∫^∞𝑡𝑎(𝑥)𝑑𝑥,

𝛽𝑏_𝑚(𝑡)=𝑎(𝑡)𝑔(𝑡,𝑚);𝑏^∗𝑚(𝑡)=𝐴(𝑡)𝑔(𝑡,𝑚),_𝑚𝛽((𝑠𝑢),=𝑠)=^∞0𝑒^−𝑠𝑡_∞𝑏_𝑚(𝑡)𝑑𝑡;𝛽^∗𝑚(𝑠)=^∞0𝑒^−𝑠𝑡𝑏^∗𝑚(𝑡)𝑑𝑡,_𝑚₌₀𝛽_𝑚(𝑠)𝑢^𝑚;𝛽^∗(𝑢,𝑠)=^∞_𝑚₌₀𝛽^∗𝑚(𝑠)𝑢^𝑚.(2.51) The Laplace transform generating function of 𝑀(𝑡) defined by 𝜋(𝑢,𝑠)=^∞0𝑒^−𝑠𝑡E𝑢^𝑀(𝑡)𝑑𝑡(2.52) is then given by 𝜋(𝑢,𝑠)=𝜑(𝑢,1,0,𝑠), so that one has from (2.50) 𝜋(𝑢,𝑠)=𝛽^∗(𝑢,𝑠)1−𝛽(𝑢,𝑠).(2.53)

This class of counting processes denoted by ACPGRP may be extended where the underlying renewal process is replaced by an MMPP or an SMMPP. We define the former as a class of age-dependent counting processes governed by an MMPP, denoted by Markov-modulated age-dependent nonhomogeneous Poisson process (MMANHPP), and the latter as a class of age-dependent counting processes governed by an SMMPP, denoted by semi-Markov-modulated age-dependent nonhomogeneous Poisson process (SMMANHPP). The two extended classes are new and become special cases of the unified counting process proposed in this paper as we will see.

All of the counting processes discussed in this section are summarized in Figure 1.

Figure 1 

Various counting processes.

3. Unified Multivariate Counting Process [𝐌(𝐭),𝐍(𝐭)]

In this section, we propose a stochastic process representing the unified multivariate counting process discussed in Section 1, which would contain all of the counting processes given in Section 2 as special cases. More specifically, we consider a system where items arrive according to an NHPP. This arrival stream is governed by a finite semi-Markov process on 𝒥={0,…,𝐽} in that the intensity function of the NHPP depends on the current state of the semi-Markov process. That is, when the semi-Markov process is in state 𝑖 with the current dwell time of 𝑥, items arrive according to a Poisson process with intensity 𝜆_𝑖(𝑥). If the semi-Markov process switches its state from 𝑖 to 𝑗, the intensity function 𝜆_𝑖(𝑥) is interrupted, the intensity function at state 𝑗 is reset to 𝜆_𝑗(0), and the arrival process resumes. Of particular interest would be the multivariate counting processes 𝑀(𝑡)^⊤=[𝑀₀(𝑡),…,𝑀_𝐽(𝑡)] and 𝑁(𝑡)=[𝑁_𝑖𝑗(𝑡)] with 𝑀_𝑖(𝑡) and 𝑁_𝑖𝑗(𝑡) counting the number of items that have arrived in state 𝑖 by time 𝑡 and the number of transitions of the semi-Markov process from state 𝑖 to state 𝑗 by time 𝑡 respectively. The two counting processes 𝑀(𝑡) and 𝑁(𝑡) enable one to introduce a variety of interesting performance indicators as we will see.

Formally, let {𝐽(𝑡)∶𝑡≥0} be a semi-Markov process on 𝒥={0,…,𝐽} governed by a matrix cumulative distribution function (c.d.f.) 𝐴(𝑥)=𝐴_𝑖𝑗(𝑥),(3.1) which is assumed to be absolutely continuous with the matrix probability density function (p.d.f.) 𝑎(𝑥)=𝑎_𝑖𝑗(𝑥)=𝑑𝑑𝑥𝐴(𝑥).(3.2) It should be noted that, if we define 𝐴_𝑖(𝑥) and 𝐴_𝑖(𝑥) by 𝐴_𝑖(𝑥)=_𝑗∈𝒥𝐴_𝑖𝑗(𝑥);𝐴_𝑖(𝑥)=1−𝐴_𝑖(𝑥),(3.3) then 𝐴_𝑖(𝑥) is an ordinary c.d.f. and 𝐴_𝑖(𝑥) is the corresponding survival function. The hazard rate functions associated with the semi-Markov process are then defined as 𝜂_𝑖𝑗(𝑥)=𝑎^𝑖𝑗(𝑥)𝐴_𝑖(𝑥),𝑖,𝑗∈𝒥.(3.4)

For notational convenience, the transition epochs of the semi-Markov process are denoted by 𝜏_𝑛,𝑛≥0, with 𝜏₀=0. The age process 𝑋(𝑡) associated with the semi-Markov process is then defined as 𝑋(𝑡)=𝑡−max𝜏_𝑛∶0≤𝜏_𝑛≤𝑡.(3.5) At time 𝑡 with 𝐽(𝑡)=𝑖 and 𝑋(𝑡)=𝑥, the intensity function of the NHPP is given by 𝜆_𝑖(𝑥). For the cumulative arrival intensity function 𝐿_𝑖(𝑥) in state 𝑖, one has 𝐿_𝑖(𝑥)=^𝑥0𝜆_𝑖(𝑦)𝑑𝑦.(3.6) The probability of observing 𝑘 arrivals in state 𝑖 within the current age of 𝑥 can then be obtained as 𝑔_𝑖(𝑥,𝑘)=𝑒^{−𝐿𝑖(𝑥)}𝐿_𝑖(𝑥)^𝑘𝑘!,𝑘=0,1,2,…,𝑖∈𝒥.(3.7)

Of interest are the multivariate counting processes 𝑀(𝑡)^⊤=𝑀₀(𝑡),…,𝑀_𝐽(𝑡),(3.8) where 𝑀_𝑖(𝑡) represents the total number of items that have arrived by time 𝑡 while the semi-Markov process stayed in state 𝑖, and 𝑁(𝑡)=𝑁_𝑖𝑗(𝑡),(3.9) with 𝑁_𝑖𝑗(𝑡) denoting the number of transitions of the semi-Markov process from state 𝑖 to state 𝑗 by time 𝑡. It is obvious that 𝑁_𝑖(𝑡)=^def∑_ℓ∈𝒥𝑁_ℓ𝑖(𝑡) denotes the number of entries into state 𝑖 by time 𝑡. The initial state is not included in 𝑁_•𝑖(𝑡) for any 𝑖∈𝒥. In other words, if 𝐽(0)=𝑖, 𝑁_•𝑖(𝑡) remains 0 until the first return of the semi-Markov process to state 𝑖. In the next section, we will analyze the dynamic behavior of [𝑀(𝑡),𝑁(𝑡)], yielding the relevant Laplace transform generating functions. In the subsequent section, all of the counting processes discussed in Section 2 would be expressed in terms of 𝑀(𝑡) and 𝑁(𝑡), thereby providing a unified approach for studying various counting processes. The associated asymptotic behavior as 𝑡→∞ would be also discussed.

4. Dynamic Analysis of [𝐌(𝐭),𝐍(𝐭)]

The purpose of this section is to examine the dynamic behavior of the multivariate stochastic process [𝑀(𝑡),𝑁(𝑡)] introduced in Section 3 by observing its probabilistic flow in the state space. Figure 2 depicts a typical sample path of the multivariate process.

Figure 2 

Typical sample path of [𝑀(𝑡),𝑁(𝑡)].

Since [𝑀(𝑡),𝑁(𝑡)] is not Markov, we employ the method of supplementary variables. More specifically, the multivariate stochastic process [𝑀(𝑡),𝑁(𝑡),𝑋(𝑡),𝐽(𝑡)] is considered. This multivariate stochastic process is Markov and has the state space 𝑆=ℤ^+𝐽+1×ℤ^{+(𝐽+1)×(𝐽+1)}×ℝ₊×𝒥, where ℤ^+𝐽+1 and ℤ^{+(𝐽+1)×(𝐽+1)} are the set of (𝐽+1) and (𝐽+1)×(𝐽+1) dimensional nonnegative integer vectors and matrices respectively, ℝ₊ is the set of nonnegative real numbers and 𝒥={0,…,𝐽}. Let 𝐹_𝑖𝑗(𝑚,𝑛,𝑥,𝑡) be the joint probability distribution of [𝑀(𝑡),𝑁(𝑡),𝑋(𝑡),𝐽(𝑡)] defined by 𝐹_𝑖𝑗𝑚,𝑛,𝑥,𝑡=P𝑀(𝑡)=𝑚,𝑁(𝑡)=𝑛,𝑋(𝑡)≤𝑥,𝐽(𝑡)=𝑗∣𝑀(0)=0,𝑁(0)=0,𝐽(0)=𝑖.(4.1)

In order to assure the differentiability of 𝐹_𝑖𝑗(𝑚,𝑛,𝑥,𝑡) with respect to 𝑥, we assume that 𝑋(0) has an absolutely continuous initial distribution function 𝐷(𝑥) with p.d.f. 𝑑(𝑥)=(𝑑/𝑑𝑥)𝐷(𝑥). (If 𝑋(0)=0 with probability 1, we consider a sequence of initial distribution functions {𝐷_𝑘(𝑥)}^∞𝑘=1 satisfying 𝐷_𝑘(𝑥)→𝑈(𝑥) as 𝑘→∞ where 𝑈(𝑥)=1 for 𝑥≥0 and 𝑈(𝑥)=0 otherwise. The desired results can be obtained through this limiting process.) One can then define 𝑓_𝑖𝑗𝑚,𝑛,𝑥,𝑡=𝜕𝜕𝑥𝐹^𝑖𝑗𝑚,𝑛,𝑥,𝑡.(4.2) By interpreting the probabilistic flow of the multivariate process [𝑀(𝑡),𝑁(𝑡),𝑋(𝑡),𝐽(𝑡)] in its state space, one can establish the following equations: 𝑓_𝑖𝑗𝑚,𝑛,𝑥,𝑡=𝛿_{𝑖=𝑗}𝛿_{{𝑚=𝑚𝑖1𝑖}}𝛿_{𝑛=0}𝑑(𝑥−𝑡)𝐴_𝑖(𝑥)𝐴_𝑖(𝑥−𝑡)𝑔^𝑖+1−𝛿_{𝑛=0}𝑡,𝑚_𝑖𝑓^𝑚𝑗_𝑘=0𝑓_𝑖𝑗𝑚−𝑘1_𝑗,𝑛,0+,𝑡−𝑥𝐴_𝑗(𝑥)𝑔_𝑗(𝑥,𝑘);(4.3)_𝑖𝑗𝑓𝑚,𝑛,0+,𝑡=1−𝛿_{𝑛=0}_ℓ∈𝒥^∞0𝑓_𝑖ℓ𝑚,𝑛−1_ℓ𝑗,𝑥,𝑡𝜂_ℓ𝑗(𝑥)𝑑𝑥;(4.4)_𝑖𝑗𝑚,𝑛,𝑥,0=𝛿_{𝑖=𝑗}𝛿_{𝑚=0}𝛿_{𝑛=0}𝑑(𝑥),(4.5) where 1_𝑖 is the column vector whose 𝑖th element is equal to 1 with all other elements being 0, 1_𝑖𝑗=1_𝑖1^⊤𝑗 and 𝑓_𝑖𝑗(𝑚,𝑛,0+,𝑡)=0for𝑁≤0.

The first term of the right-hand side of (4.3) represents the case that 𝐽(𝑡) has not left the initial state 𝐽(0)=𝑖 by time 𝑡[𝛿_{𝑖=𝑗}=1and𝛿_{𝑛=0}=1] and there have been 𝑚_𝑖 items arrived during time 𝑡[𝛿_{{𝑚=𝑚𝑖1𝑖}}=1], provided that 𝐽(0)=𝑖 and 𝑋(0)=𝑥−𝑡. The second term corresponds to the case that 𝐽(𝑡) made at least one transition from 𝐽(0)=𝑖 by time 𝑡[𝛿_{𝑛=0}=0], the multivariate process [𝑀(𝑡),𝑁(𝑡),𝑋(𝑡),𝐽(𝑡)] just entered the state [𝑚−𝑘1_𝑗,𝑛,0+,𝑗] at time 𝑡−𝑥, 𝐽(𝑡) remained in state 𝑗 until time 𝑡 with 𝑋(𝑡)=𝑥, and there have been 𝑘 items arrived during the current age 𝑥, provided that 𝐽(0)=𝑖 and 𝑋(0)=𝑥−𝑡. For the multivariate process at [𝑚,𝑛,0+,𝑗] at time 𝑡, it must be at [𝑚,𝑛−1_ℓ𝑗,𝑥,ℓ] followed by a transition from ℓ to 𝑗 at time 𝑡 which increases 𝑁_ℓ𝑗(t) by one, explaining (4.4). Equations (4.5) merely describe the initial condition that [𝑀(0),𝑁(0),𝑋(0),𝐽(0)]=[0,0,𝑥,𝑖].

In what follows, the dynamic behavior of the multivariate process [𝑀(𝑡),𝑁(𝑡),𝑋(𝑡),𝐽(𝑡)] would be captured by establishing the associated Laplace transform generating functions based on (4.3), (4.4) and (4.5). For notational convenience, the following matrix functions are employed: b𝑏_𝑘(𝑡)=𝑏_{𝑘∶𝑖𝑗}(𝑡);𝑏_{𝑘∶𝑖𝑗}(𝑡)=𝑎_𝑖𝑗(𝑡)𝑔_𝑖(𝑡,𝑘),(4.6)_{∗𝑘∶𝐷}(𝑡)=𝛿_{𝑖=𝑗}𝑏^{∗𝑘∶𝑖}(𝑡)=⎡⎢⎢⎣𝑏^∗𝑘∶0(𝑡)⋱𝑏_{∗𝑘∶𝐽}(𝑡)⎤⎥⎥⎦;𝑏^{∗𝑘∶𝑖}𝑟(𝑡)=𝐴_𝑖(𝑡)𝑔_𝑖(𝑡,𝑘),(4.7)_∗𝑘(𝑡)=𝑟^{∗𝑘∶𝑖𝑗}(𝑡);𝑟^{∗𝑘∶𝑖𝑗}(𝑡)=𝑔_𝑖(𝑡,𝑘)^∞0𝑑(𝑥−𝑡)𝑎^𝑖𝑗(𝑥)𝛽𝐴_𝑖(𝑥−𝑡)𝑑𝑥,(4.8)_𝑘𝛽(𝑠)=𝛽_{𝑘∶𝑖𝑗}(𝑠);𝛽_{𝑘∶𝑖𝑗}(𝑠)=^∞0𝑒^−𝑠𝑡𝑏_{𝑘∶𝑖𝑗}(𝑡)𝑑𝑡,(4.9)𝛽𝑢,𝑠=𝛽_𝑖𝑗𝑢_𝑖,𝑠;𝛽_𝑖𝑗𝑢_𝑖,𝑠=^∞_𝑘=0𝛽_{𝑘∶𝑖𝑗}(𝑠)𝑢^𝑘𝑖,(4.10)_{∗𝑘∶𝐷}(𝑠)=⎡⎢⎢⎣𝛽^∗𝑘∶0(𝑠)⋱𝛽_{∗𝑘∶𝐽}(𝑠)⎤⎥⎥⎦;𝛽^{∗𝑘∶𝑖}𝛽(𝑠)=^∞0𝑒^−𝑠𝑡𝑏^{∗𝑘∶𝑖}(𝑡)𝑑𝑡,(4.11)_∗𝐷𝑢,𝑠=⎡⎢⎢⎣𝛽^∗0𝑢₀,𝑠⋱𝛽_∗𝐽𝑢_𝐽,𝑠⎤⎥⎥⎦;𝛽^∗𝑖𝜌𝑢_𝑖,𝑠=^∞_𝑘=0𝛽^{∗𝑘∶𝑖}(𝑠)𝑢^𝑘𝑖,(4.12)_∗𝑘𝜌(𝑠)=𝜌^{∗𝑘∶𝑖𝑗}(𝑠);𝜌^{∗𝑘∶𝑖𝑗}(𝑠)=^∞0𝑒^−𝑠𝑡𝑟^{∗𝑘∶𝑖𝑗}(𝑡)𝑑𝑡,(4.13)_∗𝑢,𝑠=𝜌^∗𝑖𝑗𝑢_𝑖,𝑠;𝜌^∗𝑖𝑗u_𝑖,𝑠=^∞_𝑘=0𝜌^{∗𝑘∶𝑖𝑗}(𝑠)𝑢^𝑘𝑖,(4.14) where 𝑢^𝑚𝑣^𝑛=∏_𝑖∈𝒥𝑢^𝑚𝑖_𝑖∏_{(𝑖,𝑗)∈𝒥×𝒥⧵{(0,0)}}𝑣^𝑛𝑖𝑗_𝑖𝑗. We are now in a position to prove the main theorem of this section.Theorem. Let 𝑋(0)=0. Then: 𝜉̂𝜑𝑢,𝑣,0+,𝑠=𝛽̃𝑢,𝑣,𝑠𝐼−𝛽̃(𝑢,𝑣,𝑠)⁻¹;(4.19)𝑢,𝑣,𝑤,𝑠=𝐼−𝛽̃(𝑢,𝑣,𝑠)⁻¹𝛽^∗𝐷𝑢,𝑤+𝑠,(4.20) where 𝛽̃(𝑢,𝑣,𝑠)=[𝑣_𝑖𝑗⋅𝛽_𝑖𝑗(𝑢_𝑖,𝑠)].

Proof. First, we assume that 𝑋(0) has a p.d.f. 𝑑(𝑥). Substituting (4.3) and (4.5) into (4.4), one sees that 𝑓_𝑖𝑗𝑚,𝑛,0+,𝑡=×1−𝛿_{{𝑛=0}ℓ∈𝒥}^∞0𝛿_{𝑖=ℓ}𝛿_{{𝑚=𝑚𝑖1𝑖}}𝛿_{{𝑛=1ℓ𝑗}}𝑑(𝑥−𝑡)𝐴_𝑖(𝑥)𝐴_𝑖(𝑥−𝑡)𝑔^𝑖+𝑡,𝑚_𝑖𝜂_ℓ𝑗(𝑥)𝑑𝑥_ℓ∈𝒥1−𝛿^{𝑛=1_ℓ𝑗^}^∞0𝑚ℓ_𝑘=0𝑓_𝑖ℓ𝑚×−𝑘1_ℓ,𝑛−1_ℓ𝑗,0+,𝑡−𝑥=1−𝛿𝐴_ℓ(𝑥_{𝑛=0})𝑔_ℓ(𝑥,𝑘)𝜂_ℓ𝑗(𝑥)𝑑𝑥𝛿^{𝑚=𝑚_𝑖¹_𝑖^}𝛿_{{𝑛=1𝑖𝑗}}𝑔_𝑖𝑡,𝑚_𝑖^∞0𝑑(𝑥−𝑡)𝑎^𝑖𝑗(𝑥)+𝐴_𝑖(𝑥−𝑡)𝑑𝑥_ℓ∈𝒥1−𝛿^{𝑛=1_ℓ𝑗^}^𝑚ℓ_𝑘=0^∞0𝑓_𝑖ℓ𝑚×𝑎_ℓ𝑗−𝑘1_ℓ,𝑛−1_ℓ𝑗,0+,𝑡−𝑥(𝑥)𝑔_ℓ(𝑥,𝑘)𝑑𝑥.(4.21) Consequently, one sees from (4.6) and (4.8) that 𝑓_𝑖𝑗𝑚,𝑛,0+,𝑡=×1−𝛿𝛿_{𝑛=0}_{𝑚=𝑚𝑖1𝑖}𝛿_{𝑛=1𝑖𝑗}𝑟^∗𝑚_{𝑖∶𝑖𝑗}(𝑡)+_ℓ∈𝒥×_𝑚ℓ1−𝛿^{𝑛=1_ℓ𝑗^}_𝑘=0^∞0𝑓_𝑖ℓ𝑚−𝑘1_ℓ,𝑛−1_ℓ𝑗,0+,𝑡−𝑥𝑏_{𝑘∶ℓ𝑗}(𝑥)𝑑𝑥.(4.22) where 𝑎_𝑖𝑗(𝑥)=𝐴_𝑖(𝑥)𝜂_𝑖𝑗(𝑥) is employed from (3.4). By taking the Laplace transform of both sides of (4.22) with respect to 𝑡, it follows that 𝜉_𝑖𝑗𝑚,𝑛,0+,𝑠=×1−𝛿𝛿_{𝑛=0}_{𝑚=𝑚𝑖1𝑖}+𝛿_{𝑛=1𝑖𝑗}𝜌^∗𝑚_{𝑖∶𝑖𝑗}(𝑠)_ℓ∈𝒥1−𝛿^{𝑛=1_ℓ𝑗^}^𝑚ℓ_𝑘=0𝜉_𝑖ℓ𝑚−𝑘1_ℓ,𝑛−1_ℓ𝑗,0+,𝑠𝛽_{𝑘∶ℓ𝑗}(𝑠).(4.23) By taking the multivariate generating function of (4.23) with respect to 𝑚 and 𝑛, it can be seen that 𝜉̂_𝑖𝑗𝑢,𝑣,0+,𝑠=_{𝑚∈ℤ+𝐽+1}_{𝑛∈ℤ+(𝐽+1)×(𝐽+1)⧵{0}}𝜉_𝑖𝑗𝑚,𝑛,0+,𝑠𝑢^𝑚𝑣^𝑛=_{𝑚∈ℤ+𝐽+1𝑛∈ℤ+(𝐽+}_{1)×(𝐽+1)⧵{0}}𝛿^{{𝑚=𝑚𝑖1𝑖}}𝛿_{{𝑛=1𝑖𝑗}}𝜌^∗𝑚_{𝑖∶𝑖𝑗}(𝑠)u^𝑚𝑣^𝑛+_ℓ∈𝒥1−𝛿^{𝑛=1_ℓ𝑗^}^𝑚ℓ_𝑘=0𝜉_𝑖ℓ𝑚−𝑘1_ℓ,𝑛−1_ℓ𝑗,0+,𝑠𝛽_{𝑘∶ℓ𝑗}(𝑠)𝑢^𝑚𝑣^𝑛.(4.24) It then follows from (4.10), (4.14), (4.16) and the discrete convolution property that 𝜉̂_𝑖𝑗𝑢,𝑣,0+,𝑠=+^∞_𝑚𝑖=0𝑢^𝑚𝑖_𝑖𝑣_𝑖𝑗𝜌^∗𝑚_{𝑖∶𝑖𝑗}(𝑠)_ℓ∈𝒥⎛⎜⎜⎜⎝𝑣_ℓ𝑗_{𝑚∈ℤ+𝐽+1}^𝑚ℓ_𝑘=0_{𝑛∈ℤ+(𝐽+1)×(𝐽+1)⧵{0}}𝜉_𝑖ℓ𝑚−𝑘1_ℓ,𝑛,0+,𝑠𝛽_{𝑘∶ℓ𝑗}(𝑠)𝑢^𝑚𝑣^𝑛=𝑣_𝑖𝑗⎞⎟⎟⎟⎠𝜌^∗𝑖𝑗𝑢_𝑖,𝑠+_ℓ∈𝒥𝑣_ℓ𝑗𝜉̂_𝑖ℓ𝑢,𝑣,0+,𝑠𝛽_ℓ𝑗𝑢_ℓ,𝑠.(4.25) The last expression can be rewritten in matrix form, and one has 𝜉̂𝑢,𝑣,0+,𝑠=𝜌̃^∗𝑢,𝑣,𝑠+𝜉̂𝑢,𝑣,0+,𝑠𝛽̃𝑢,𝑣,𝑠,(4.26) which can be solved for 𝜉̂(𝑢,𝑣,0+,𝑠) as 𝜉̂𝑢,𝑣,0+,𝑠=𝜌̃^∗𝑢,𝑣,𝑠𝐼−𝛽̃(𝑢,𝑣,𝑠)⁻¹,(4.27) where 𝜌̃^∗(𝑢,𝑣,s)=[𝑣_𝑖𝑗⋅𝜌^∗𝑖𝑗(𝑢_𝑖,𝑠)].
Next, we introduce the following double Laplace transform:
𝜀_𝑘∶𝑖(𝑤,𝑠)=^∞0𝑒^−𝑤𝑥𝑒^−𝑠𝑡𝑑(𝑥−𝑡)𝐴_𝑖(𝑥)𝐴_𝑖(𝑥−𝑡)𝑔^𝑖(𝑡,𝑘)𝑑𝑡𝑑𝑥,(4.28) and the associated diagonal matrix 𝜀_𝐷𝑢,𝑤,𝑠=⎡⎢⎢⎣𝜀^∗0𝑢₀,𝑤,𝑠⋱𝜀_∗𝐽𝜀𝑢_𝐽,𝑤,𝑠_𝑖⎤⎥⎥⎦;𝑢_𝑖,𝑤,𝑠=^∞_𝑘=0𝜀_𝑘∶𝑖(𝑤,𝑠)𝑢^𝑘𝑖.(4.29) By taking the double Laplace transform of (4.3), one sees from (4.7) and (4.28) that 𝜑_𝑖𝑗𝑚,𝑛,𝑤,𝑠=+𝛿1_{𝑖=𝑗}−𝛿𝛿_{{𝑚=𝑚𝑖1𝑖}{𝑛=0}}𝛿_{𝑛=0}𝜀_{𝑚𝑖∶𝑖}(𝑤,𝑠)^𝑚𝑗_𝑘=0𝜉_𝑖𝑗𝑚−𝑘1_𝑗,𝑛,0+,𝑠𝛽^{∗𝑘∶𝑗}(𝑤+𝑠).(4.30) By taking the double generating function, this then leads to 𝜑_𝑖𝑗=𝑢,𝑣,𝑤,𝑠_{𝑚∈ℤ+𝐽+1}_𝑛∈ℤ+(_{𝐽+1)×(𝐽+1)}𝜑_𝑖𝑗𝑚,𝑛,𝑤,𝑠𝑢^𝑚𝑣^𝑛=𝛿_{∞{𝑖=𝑗}𝑚}_𝑖=0+𝑢^𝑚𝑖_𝑖𝜀_{𝑚𝑖∶𝑖}(𝑤,𝑠)_{𝑚∈ℤ+𝐽+1𝑛∈ℤ+(𝐽+1}_{)×(𝐽+1)⧵0}^𝑚𝑗_𝑘=0𝜉_𝑖𝑗𝑚−𝑘1_𝑗,𝑛,0+,𝑠𝛽^{∗𝑘∶𝑗}(𝑤+𝑠)𝑢^𝑚𝑣^𝑛=𝛿_{𝑖=𝑗}𝜀_𝑖𝑢_𝑖,𝑤,𝑠+𝜉̂_𝑖𝑗𝑢,𝑣,0+,𝑠𝛽^∗𝑗𝑢_𝑗,𝑤+𝑠.(4.31) By rewriting the last expression in matrix form, it can be seen that 𝜑𝑢,𝑣,𝑤,𝑠=𝜀_𝐷𝑢,𝑤,𝑠+𝜉̂𝑢,𝑣,0+,𝑠𝛽^∗𝐷𝑢,𝑤+𝑠.(4.32)
We now consider the limiting process 𝐷(𝑥)→𝑈(𝑥). For the p.d.f., this limiting process becomes 𝑑(𝑥)→𝛿(𝑥) where 𝛿(𝑥) is the Delta function defined as a unit function for convolution, that is, 𝑓(𝑥)=∫𝑓(𝑥−𝜏)𝛿(𝜏)𝑑𝜏 for any integrable function 𝑓. Accordingly, it can be seen from (4.8) that 𝑟^{∗𝑘∶𝑖𝑗}(𝑡)→𝑏_{𝑘∶𝑖𝑗}(𝑡). This in turn implies from (4.28) that 𝜀_𝑘∶𝑖(𝑤,𝑠)→𝛽^{∗𝑘∶𝑖}(𝑤+𝑠). Consequently, it follows in matrix form that 𝜌̃^∗(𝑢,𝑣,𝑠)→𝛽̃(𝑢,𝑣,𝑠) and 𝜀_𝐷(𝑢,𝑤,𝑠)→𝛽^∗𝐷(𝑢,𝑤+𝑠). From (4.27), it can be readily seen that 𝜉̂(𝑢,𝑣,0+,𝑠)→𝛽̃(𝑢,𝑣,𝑠){𝐼−𝛽̃(𝑢,𝑣,𝑠)}⁻¹, proving (4.19). One also sees from (4.28) that
𝜑𝑢,𝑣,𝑤,𝑠⟶=𝐼𝛽^∗𝐷𝑢,𝑤+𝑠+𝜉̂𝑢,𝑣,0+,𝑠𝛽^∗𝐷𝑢,𝑤+𝑠+𝛽̃𝑢,𝑣,𝑠𝐼−𝛽̃(𝑢,𝑣,𝑠)⁻¹𝛽^∗𝐷=𝐼𝑢,𝑤+𝑠=𝐼−𝛽̃𝑢,𝑣,𝑠+𝛽̃𝑢,𝑣,𝑠𝐼−𝛽̃(𝑢,𝑣,𝑠)⁻¹𝛽^∗𝐷𝑢,𝑤+𝑠−𝛽̃(𝑢,𝑣,𝑠)⁻¹𝛽^∗𝐷𝑢,𝑤+𝑠,(4.33) which proves (4.20), completing the proof.In the next section, it will be shown that all the transform results obtained in Section 2 can be derived from Theorem 4.1.

5. Derivation of the Special Cases from the Unified Counting Process

We are now in a position to demonstrate the fact that the proposed multivariate counting process introduced in Section 3 and analysed in Section 4 indeed unifies the existing counting processes discussed in Section 2. We will do so by deriving the transform results of Section 2 from Theorem 4.1.

5.1. Derivation of Poisson Process

Let 𝑁(𝑡) be a Poisson process with intensity 𝜆 as discussed in Section 2.1. From (2.5), one sees that 𝜋(𝑣,𝑠)=∫^∞0𝑒^−𝑠𝑡E[𝑣^𝑁(𝑡)]𝑑𝑡=∫^∞0𝑒^−𝑠𝑡𝜋(𝑣,𝑡)𝑑𝑡 is given by 𝜋(𝑣,𝑠)=1𝑠+𝜆(1−𝑣).(5.1)

For the unified counting process, we consider a single state Markov-chain in continuous time where only the number of self transitions in (0,𝑡] is counted. More specifically, let 𝒥={0}, 𝑁(𝑡)=[𝑁₀₀(𝑡)], 𝑢₀=1, 𝑣₀₀=𝑣, 𝑤=0, 𝜆₀(𝑡)=0 and 𝑎₀₀(𝑥)=𝜆𝑒^−𝜆𝑥. We note from (3.7) that 𝜆₀(𝑡)=0 implies 𝑔₀(𝑥,𝑘)=𝛿_{𝑘=0} so that 𝑏_𝑘∶00(𝑡)=𝛿_{𝑘=0}𝑎₀₀(𝑡) from (4.6). This then implies 𝛽₀₀(1,𝑠)=𝜆/(𝑠+𝜆). Similarly, since 𝐴₀(𝑥)=𝑒^−𝜆𝑥, one has 𝛽^∗0(1,𝑠)=1/(𝑠+𝜆). It then follows from Theorem 4.1 that 𝜑₀(1,𝑣,0+,𝑠)=11−𝑣𝛽₀₀(1,𝑠)𝛽^∗0(1,𝑠)=1𝑠+𝜆(1−𝑣).(5.2) Hence, from (5.1) and (5.2), one concludes that 𝜋(𝑣,𝑠)=𝜑₀(1,𝑣,0+,𝑠).

5.2. Derivation of NHPP

Let 𝑀(𝑡) be an NHPP of Section 2.2 characterized by a time dependent intensity function 𝜆(𝑡). It can be seen from (2.8) that 𝜋(𝑢,𝑠)=∫^∞0𝑒^−𝑠𝑡E[𝑢^𝑀(𝑡)]𝑑𝑡=∫^∞0𝑒^−𝑠𝑡𝜋(𝑢,𝑡)𝑑𝑡 is given by 𝜋(𝑢,𝑠)=^∞0𝑒^−𝑠𝑡𝑒^{−𝐿(𝑡)(1−𝑢)}𝑑𝑡.(5.3)

In order to see that 𝑀(𝑡) can be viewed as a special case of the proposed multivariate counting process, we first consider a single state semi-Markov process where the dwell time in the state is deterministic given by 𝑇. The marginal counting process 𝑀(𝑡)=[𝑀₀(𝑡)] then converges in distribution to 𝑀(𝑡) as 𝑇→∞ as we show next.

Let 𝒥={0}, 𝑢₀=𝑢, 𝑣₀₀=1, 𝑤=0, and 𝜆₀(𝑥)=𝜆(𝑥). It then follows that 𝑎₀₀(𝑥)=𝛿(𝑥−𝑡) and therefore 𝑏_𝑘∶00(𝑡)=𝛿(𝑡−𝑇)𝑔₀(𝑡,𝑘) from (4.6). This in turn implies that 𝛽₀₀(𝑢,𝑠)==𝑒^∞0𝑒^−𝑠𝑡𝛿(𝑡−𝑇_{−{𝑠𝑇+𝐿0(𝑇)(1−𝑢)}})𝑒^{−𝐿0(𝑡)(1−𝑢)}𝑑𝑡.(5.4) Let 𝑈(𝑡) be the step function defined by 𝑈(𝑡)=0 if 𝑡<0 and 𝑈(𝑡)=1 otherwise. Since 𝐴₀₀(𝑡)=1−𝑈(𝑡−𝑇)=𝛿_{{0≤𝑡<𝑇}}, one sees that 𝛽^∗0(𝑢,𝑠)==^∞0𝑒^−𝑠𝑡𝛿_{{0≤𝑡<𝑇}}𝑒^{−𝐿0(𝑡)(1−𝑢)}𝑑𝑡_𝑇0𝑒^−𝑠𝑡𝑒^{−𝐿0(𝑡)(1−𝑢)}𝑑𝑡.(5.5) Theorem 4.1 together with (5.4) and (5.5) then leads to 𝜑₀₀(𝑢,1,0+,𝑠)=11−𝛽₀₀(𝑢,𝑠)𝛽^∗0(𝑢,𝑠)=∫^𝑇0𝑒^−𝑠𝑡𝑒^{−𝐿0(𝑡)(1−𝑢)}𝑑𝑡1−𝑒^{−{𝑠𝑇+𝐿0(𝑇)(1−𝑢)}}.(5.6) Now it can be readily seen that 𝜑₀₀(𝑢,1,0+,𝑠) in (5.6) converges to 𝜋(𝑢,𝑠) in (5.3) as 𝑇→∞, proving the claim.

5.3. Derivation of MMPP

Let 𝑀(𝑡) be an MMPP of Section 2.3 characterized by (𝜈,𝜆_𝐷). We show that the Laplace transform generating function 𝜋(𝑢,𝑠)=∫^∞0𝑒^−𝑠𝑡E[𝑢^𝑀(𝑡)]𝑑𝑡 given in (2.16) can be derived as a special case of Theorem 4.1.

For the unified multivariate counting process, let 𝒥={0,…,𝐽}, 𝑀(𝑡)=∑_𝑖∈𝒥𝑀_𝑖(𝑡), 𝑢=𝑢1, 𝑣=1, 𝑤=0, 𝜆_𝑖(𝑡)=𝜆_𝑖, and 𝑎(𝑥)=[𝑒^{−𝜈𝑖𝑥}𝜈_𝑖𝑗]=𝑒^{−𝜈𝐷𝑥}𝜈. From (3.3) one sees that 𝐴_𝐷(𝑥)=𝐼−𝐴_𝐷=[𝛿_{𝑖=𝑗}{1−∑_𝑗∈𝒥𝐴_𝑖𝑗(𝑥)}]=[𝛿_{𝑖=𝑗}𝑒^{−𝜈𝑖𝑥}]=𝑒^{−𝜈𝐷𝑥}. Therefore, one sees from (4.6), (4.9) and (4.10) that 𝛽𝑢1,𝑠=^∞_𝑘=0^∞0𝑒^−𝑠𝑡𝑔_𝐷(𝑡,𝑘)𝑎(𝑡)𝑑𝑡𝑢^𝑘1=𝑠𝐼+𝜈_𝐷+(1−𝑢)𝜆_𝐷⁻¹𝜈(5.7) and similarly one has, from (4.7), (4.11) and (4.12), 𝛽^∗𝐷𝑢1,𝑠=^∞_𝑘=0^∞0𝑒^−𝑠𝑡𝐴_𝐷(𝑡)𝑔_𝐷(𝑡,𝑘)𝑑𝑡𝑢^𝑘1=𝑠𝐼+𝜈_𝐷+(1−𝑢)𝜆_𝐷⁻¹.(5.8) It then follows from Theorem 4.1, (5.7)and (5.8) that 𝜑𝑢1,1,0+,𝑠==𝐼𝐼−𝛽̃(𝑢1,1,𝑠)⁻¹𝛽^∗𝐷𝑢1,𝑠=𝑠𝐼−𝛽(𝑢1,𝑠)⁻¹𝛽^∗𝐷𝑢1,𝑠+𝜈_𝐷+(1−𝑢)𝜆_𝐷−𝜈⁻¹=𝑠𝐼−𝑄+(1−𝑢)𝜆_𝐷⁻¹,(5.9) which coincides with (2.16) as expected.

5.4. Derivation of Renewal Process

In order to demonstrate that a renewal process is a special case of the unified multivariate counting process, we follow the line of the arguments for the case of Poisson processes. Namely, let 𝒥={0}, 𝑁(𝑡)=𝑁₀(𝑡), 𝑢₀=1, 𝑣₀₀=𝑣, 𝑤=0, 𝑎₀(𝑥)=𝑎(𝑥) and 𝐴₀(𝑥)=1−∫^𝑥0𝑎(𝑦)𝑑𝑦. From Theorem 4.1, one has 𝜑₀(1,𝑣,0+,𝑠)=11−𝑣𝛽(1,𝑠)𝛽^∗=1(1,𝑠)1−𝑣∫^∞0𝑒^−𝑠𝑡𝑎(𝑡)𝑑𝑡×^∞0=1𝑒^−𝑠𝑡𝐴(𝑡)𝑑𝑡1−𝑣𝛼(𝑠)⋅1−𝑠𝛼(𝑠),(5.10) which agrees with (2.19).

5.5. Derivation of MRP

Let 𝑁_ℓ(𝑡)=[𝑁_ℓ0(𝑡),…,𝑁_ℓ𝐽(𝑡)] be an MRP discussed in Section 2.5. We recall that 𝑁_ℓ𝑟(𝑡) describes the number of entries of the semi-Markov process into state 𝑟 in (0,𝑡] given that 𝐽(0)=ℓ. For the unified multivariate counting process, 𝑁_𝑖𝑗(𝑡) counts the number of transitions from state 𝑖 to state 𝑗 in (0,𝑡]. Hence, one has 𝑁_ℓ𝑟(𝑡)=∑_𝑖∈𝒥𝑁_𝑖𝑟(𝑡) provided that 𝐽(0)=ℓ. Accordingly, we set 𝑣=[𝑣₀̃1,…,𝑣_𝑟̃1,…,𝑣_𝐽̃1], that is, 𝑣_𝑖𝑟=𝑣_𝑟̃ for all 𝑖∈𝒥. With 𝑤=0+, 𝑢=1, 𝜆_ℓ(𝑡)=0 for all ℓ∈𝒥, one has 𝛽_ℓ𝑟(1,𝑠)=𝛼_ℓ𝑟(𝑠) and 𝛽^∗ℓ={1−𝛼_ℓ(𝑠)}/𝑠 from (4.6), (4.7) and (4.9) through (4.12), where 𝛼_ℓ(𝑠)=∑_𝑟∈𝒥𝛼_ℓ𝑟(𝑠). Let 𝑣_𝐷̃=[𝛿_{ℓ=𝑟}𝑣_𝑟̃]. It then follows from Theorem 4.1 that 𝜑1,𝑣,0+,𝑠=𝐼−𝑣_r̃𝛼_ℓ𝑟(𝑠)⁻¹×𝛿_{ℓ=𝑟}1−𝛼_ℓ(𝑠)=𝐼𝑠−𝛼(𝑠)𝑣_𝐷̃⁻¹×𝛿_{ℓ=𝑟}1−𝛼_ℓ(𝑠)𝑠.(5.11) It should be noted that, with 𝑣̃=[𝑣₀̃,…,𝑣_𝐽̃]^⊤ and 𝑣̃𝑁_(𝑡)=∏_𝑟∈𝒥𝑣𝑁̃^{ℓ𝑟𝑟(𝑡)}, the ℓ−𝑟 element of 𝜑(1,𝑣,0+,𝑠) in (5.11) can be written as 𝜑_ℓ𝑟1,𝑣,0+,𝑠=ℒE𝑣̃𝑁_(𝑡),𝐽(𝑡)=𝑟∣𝐽(0)=ℓ.(5.12)

We now focus on 𝑁_ℓ𝑟(𝑡) for ℓ=0,1,…,𝐽. In doing so, let 𝑣̃=_{𝐷∶𝑗def}[1₀,…,𝑣_𝑗̃1_𝑗,…,1_𝐽] and define 𝜓𝑣_𝑟̃,𝑠=^def⎡⎢⎢⎢⎢⎢⎣ℒE𝑣𝑁̃^{𝑟𝑟(𝑡)}ℒE𝑣𝑁̃⋮∣𝐽(0)=0_{𝑟𝑟(𝑡)}∣𝐽(0)=𝐽⎤⎥⎥⎥⎥⎥⎦.(5.13) It then follows from (5.11) through (5.13) that 𝜓𝑣_𝑟̃,𝑠=𝜑1,𝑣̃_𝐷∶𝑟,0+,𝑠1=𝐼−𝛼(𝑠)𝑣̃_𝐷∶𝑟⁻¹×𝛿_{ℓ=𝑟}1−𝛼_ℓ(𝑠)𝑠×1,(5.14) that is, one has 𝜓𝑣_𝑟̃,𝑠=1𝑠𝐼−𝛼(𝑠)𝑣̃_𝐷∶𝑟⁻¹×1−𝛼(𝑠)1.(5.15)

We recall that 𝐻_ℓ𝑟(𝑡)=E[𝑁_ℓ𝑟(𝑡)], which can be obtained by differencing 𝜓(𝑣_𝑟̃,𝑠) with respect to 𝑣_𝑟̃ at 𝑣_𝑟̃=1. More formally, one has 𝜕𝜕𝑣_𝑟̃𝜓(𝑣_𝑟̃,𝑠)|||𝑣̃_𝑟=1=⎡⎢⎢⎢⎣ℒE𝑁^𝑟ℒE𝑁(𝑡)⋮∣𝐽(0)=0_𝑟(𝑡)∣𝐽(0)=𝐽⎤⎥⎥⎥⎦,(5.16) which is the 𝑟th column of ℒ{𝐻(𝑡)} given in (2.20). By noting that 𝑑𝑑𝑥𝐼−𝑓(𝑥)⁻¹=𝐼−𝑓(𝑥)⁻¹𝑑𝑑𝑥𝑓(𝑥)𝐼−𝑓(𝑥)⁻¹,(5.17) one sees from (5.15) that 𝜕𝜕𝑣_𝑗̃𝜓(𝑣_𝑗̃,𝑠)|||𝑣̃_𝑗=1=1𝑠𝐼−𝛼(𝑠)⁻¹⎡⎢⎢⎣0𝛼^0𝑗(𝑠)⋮0𝛼_𝐽𝑗(𝑠)⎤⎥⎥⎦𝐼−𝛼(𝑠)⁻¹1−𝛼(𝑠)1.(5.18) Since {𝐼−𝛼(𝑠)}⁻¹=∑^∞𝑘=0𝛼(𝑠)^𝑘 and 𝛼⁰(𝑠)=𝐼, it can be seen that 𝐼−𝛼(𝑠)⁻¹1−𝛼(𝑠)1=^∞_𝑘=0𝛼(𝑠)^𝑘1−𝛼(𝑠)1=1.(5.19) Substituting (5.19) into (5.18), one then concludes that 𝜕𝜕𝑣_𝑗̃𝜓(𝑣_𝑗̃,𝑠)|||𝑣̃_𝑗=1=1𝑠𝐼−𝛼(𝑠)⁻¹⎡⎢⎢⎣𝛼𝛼⋮_0𝑗(𝑠)_𝐽𝑗(𝑠)⎤⎥⎥⎦.(5.20) This in turn implies that 𝜕𝜕𝑣₀̃𝜓(𝑣₀̃,𝑠)||||𝑣̃₀₌₁,…,𝜕𝜕𝑣_𝐽̃𝜓(𝑣_𝐽̃,𝑠)||||𝑣̃_𝐽=1=1𝑠𝐼−𝛼(𝑠)⁻¹𝛼(𝑠),(5.21) which agrees with ℒ{𝐻(𝑡)} of (2.20), completing the derivation.

5.6. Derivation of NESMPS

As in Section 2.6, let the state space 𝒥 of 𝐽(𝑡) be decomposed into a set of good states 𝐺(≠𝜙) and a set of bad states 𝐵(≠𝜙) satisfying 𝐽=𝐺∪𝐵 and 𝐺∩𝐵=𝜙. The counting process 𝑁_𝐺𝐵(𝑡) of Section 2.6 describes the number of entries of 𝐽(𝑡) into 𝐵 by time 𝑡. The Laplace transform generating function of the joint probability of 𝑁_𝐺𝐵(𝑡), the age process 𝑋(𝑡) and 𝐽(𝑡) is given in (2.26).

In the context of the unified multivariate counting process [𝑀(𝑡),𝑁(𝑡)] discussed in Section 3, one expects to have 𝑁_𝐺𝐵(𝑡)=_𝑖∈𝐺_𝑗∈𝐵𝑁_𝑖𝑗(𝑡).(5.22) In order to prove (5.22) formally, we set 𝑢=1;𝑣=1^𝐵𝐵1_𝐵𝐺𝑣1_𝐺𝐵1_𝐺𝐺.(5.23) From (3.7), (4.6), (4.9) and (4.10), one has 𝛽_𝑖𝑗(1,𝑠)=𝛼_𝑖𝑗(𝑠) so that 𝛽̃1,𝑣,𝑠=𝛼^𝐵𝐵𝑣𝛼_𝐺𝐵(𝑠)𝛼_𝐵𝐺(𝑠)(𝑠)𝛼_𝐺𝐺(𝑠),(5.24) where 𝛽̃(𝑢,𝑣,𝑠) is as given in Theorem 4.1. Similarly, it can be seen from (3.7), (4.7), (4.11) and (4.12) that 𝛽^∗𝑖(1,𝑤+𝑠)={1−𝛼_𝑖(𝑤+𝑠)}/(𝑤+𝑠) and hence 𝛽^∗𝐷1,𝑤+𝑠=1𝑤+𝑠𝐼−𝛼_𝐷(𝑤+𝑠).(5.25) Substituting (5.24) and (5.25) into (4.20), it then follows that 𝜑1,𝑣,𝑤,𝑠=1𝑤+𝑠𝐼−𝛽̃(1,𝑣,𝑠)⁻¹𝐼−𝛼_𝐷(𝑤+𝑠).(5.26)

By comparing (2.26) with (5.26), (5.22) holds true if and only if 𝛾₀(𝑠)𝐼−𝑣𝛽(𝑠)⁻¹=𝐼−𝛽̃(1,𝑣,𝑠)⁻¹.(5.27) In what follows, we prove that (5.27) indeed holds true. From (5.24), one sees that 𝐼−𝛽̃1,𝑣,𝑠=⎡⎢⎢⎣𝜒^𝐵−1−𝑣𝛼_𝐺𝐵(𝑠)−𝛼_𝐵𝐺(𝑠)(𝑠)𝜒^𝐺−1(𝑠)⎤⎥⎥⎦,(5.28) where 𝜒_𝐺(𝑠) and 𝜒_𝐵(𝑠) are as given in (2.28). We define the inverse matrix of (5.28) by 𝐼−𝛽̃(1,𝑣,𝑠)⁻¹=𝐶^𝐵𝐵𝐶(𝑣,𝑠)𝐶_𝐵𝐺(𝑣,𝑠)_𝐺𝐵(𝑣,𝑠)𝐶_𝐺𝐺(𝑣,𝑠).(5.29) Since the multiplication of the two matrices in (5.28) and (5.29) yields the identity matrix, it follows that

𝜒^𝐵−1(𝑠)𝐶_𝐵𝐵(𝑣,𝑠)−𝛼_𝐵𝐺(𝑠)𝐶_𝐺𝐵(𝑣,𝑠)=𝐼_𝐵𝐵𝜒^𝐵−1(𝑠)𝐶_𝐵𝐺(𝑣,𝑠)−𝛼_𝐵𝐺(𝑠)𝐶_𝐺𝐺(𝑣,𝑠)=𝐼_𝐵𝐺−𝑣𝛼_𝐺𝐵(𝑠)𝐶_𝐵𝐵(𝑣,𝑠)+𝜒^𝐺−1(𝑠)𝐶_𝐺𝐵(𝑣,𝑠)=𝐼_𝐺𝐵−𝑣𝛼_𝐺𝐵(𝑠)𝐶_𝐵𝐺(𝑣,𝑠)+𝜒^𝐺−1(𝑠)𝐶_𝐺𝐺(𝑣,𝑠)=𝐼_𝐺𝐺.(5.30) Solving the above equations for 𝐶_••(𝑣,𝑠), one has

𝐶_𝐵𝐵(𝑣,𝑠)=×𝜒_𝐵(𝐼𝑠)+𝑣𝜒_𝐵(𝑠)𝛼_𝐵𝐺(𝑠)𝜒_𝐺(𝑠)_𝐺𝐺𝐶−𝑣𝛼_𝐺𝐵(𝑠)𝜒_𝐵(𝑠)𝛼_𝐵𝐺(𝑠)𝜒_𝐺(𝑠)⁻¹𝛼_𝐺𝐵(𝑠)𝜒_𝐵(𝑠),(5.31)_𝐵𝐺(𝑣,𝑠)=×𝜒_𝐵(𝐼𝑠)𝛼_𝐵𝐺(𝑠)𝜒_𝐺(𝑠)_𝐺𝐺𝐶−𝑣𝛼_𝐺𝐵(𝑠)𝜒_𝐵(𝑠)𝛼_𝐵𝐺(𝑠)𝜒_𝐺(𝑠)⁻¹,(5.32)_𝐺𝐵(𝑣,𝑠)=×𝑣𝜒_𝐺𝐼(𝑠)_𝐺𝐺𝐶−𝑣𝛼_𝐺𝐵(𝑠)𝜒_𝐵(𝑠)𝛼_𝐵𝐺(𝑠)𝜒_𝐺(𝑠)⁻¹𝛼_𝐺𝐵(𝑠)𝜒_𝐵(𝑠),(5.33)_𝐺𝐺(𝑣,𝑠)=𝜒_𝐺(𝑠)𝐼_𝐺𝐺−𝑣𝛼_𝐺𝐵(𝑠)𝜒_𝐵(𝑠)𝛼_𝐵𝐺(𝑠)𝜒_𝐺(𝑠)⁻¹.(5.34)

We next turn our attention to the left-hand side of (5.27). From (2.29), one sees that 𝐼−𝑣𝛽(𝑠)=𝐼^𝐵𝐵0_𝐵𝐺−𝑣𝛼_𝐺𝐵(𝑠)𝜒_𝐵(𝑠)𝐼_𝐺𝐺−𝑣𝛼_𝐺𝐵(𝑠)𝜒_𝐵(𝑠)𝛼_𝐵𝐺(𝑠)𝜒_𝐺(𝑠).(5.35) As before, we define the inverse matrix of (5.35) as 𝐼−𝑣𝛽(𝑠)⁻¹=𝐷^𝐵𝐵𝐷(𝑣,𝑠)𝐷_𝐵𝐺(𝑣,𝑠)_𝐺𝐵(𝑣,𝑠)𝐷_𝐺𝐺(𝑣,𝑠).(5.36) Multiplying the two matrices in (5.35) and (5.36) then yields 𝐷_𝐵𝐵(𝑣,𝑠)=𝐼_𝐵𝐵𝐷_𝐵𝐺(𝑣,𝑠)=0_𝐵𝐺−𝑣𝛼_𝐺𝐵(𝑠)𝜒_𝐵(𝑠)𝐷_𝐵𝐵(𝑣,𝑠)+𝐼_𝐺𝐺−𝑣𝛼_𝐺𝐵(𝑠)𝜒_𝐵(𝑠)𝛼_𝐵𝐺(𝑠)𝜒_𝐺(𝑠)𝐷_𝐺𝐵(𝑣,𝑠)=0_𝐺𝐵−𝑣𝛼_𝐺𝐵(𝑠)𝜒_𝐵(𝑠)𝐷_𝐵𝐺(𝑣,𝑠)+𝐼_𝐺𝐺−𝑣𝛼_𝐺𝐵(𝑠)𝜒_𝐵(𝑠)𝛼_𝐵𝐺(𝑠)𝜒_𝐺(𝑠)𝐷_𝐺𝐺(𝑣,𝑠)=𝐼_𝐺𝐺,(5.37) which in turn can be solved for 𝐷_••(𝑣,𝑠) as

𝐷𝐷_𝐵𝐵(𝑣,𝑠)=𝐼_𝐵𝐵,(5.38)_𝐵𝐺𝐷(𝑣,𝑠)=0_𝐵𝐺,(5.39)_𝐺𝐵(𝑣𝐷,𝑠)=𝑣𝐼_𝐺𝐺−𝑣𝛼_𝐺𝐵(𝑠)𝜒_𝐵(𝑠)𝛼_𝐵𝐺(𝑠)𝜒_𝐺(𝑠)⁻¹𝛼_𝐺𝐵(𝑠)𝜒_𝐵(𝑠),(5.40)_𝐺𝐺(𝑣,𝑠)=𝐼_𝐺𝐺−𝑣𝛼_𝐺𝐵(𝑠)𝜒_𝐵(𝑠)𝛼_𝐵𝐺(𝑠)𝜒_𝐺(𝑠)⁻¹.(5.41) Let the left-hand side matrix of (5.27) be described as 𝛾₀(𝑠)𝐼−𝑣𝛽(𝑠)⁻¹=𝑍^𝐵𝐵𝑍(𝑣,𝑠)𝑍_𝐵𝐺(𝑣,𝑠)_𝐺𝐵(𝑣,𝑠)𝑍_𝐺𝐺(𝑣,𝑠).(5.42) From (2.30) and (5.36) through (5.41), one sees that 𝑍_𝐵𝐵(𝑣,𝑠)=×𝐼_𝐵𝐵𝐼+𝛼_𝐵𝐵(𝑠)𝜒_𝐵(𝑠)+𝑣𝜒_𝐵(𝑠)𝛼_𝐵𝐺(𝑠)𝜒_𝐺(𝑠)_𝐺𝐺−𝑣𝛼_𝐺𝐵(𝑠)𝜒_𝐵(𝑠)𝛼_𝐵𝐺(𝑠)𝜒_𝐺(𝑠)⁻¹𝛼_𝐺𝐵(𝑠)𝜒_𝐵(𝑠).(5.43) From (2.28), one easily sees that 𝜒_𝐵(𝑠)=𝐼_𝐵𝐵+𝛼_𝐵𝐵(𝑠)𝜒_𝐵(𝑠), and hence 𝑍_𝐵𝐵(𝑣,𝑠)=𝐶_𝐵𝐵(𝑣,𝑠) from (5.31). The fact that 𝑍_𝐵𝐺(𝑣,𝑠)=𝐶_𝐵𝐺(𝑣,𝑠) is straightforward from (5.32). With 𝜒_𝐺(𝑠)=𝐼_𝐺𝐺+𝛼_𝐺𝐺(𝑠)𝜒_𝐺(𝑠), one has 𝑍_𝐺𝐵(𝑣,𝑠)=𝐶_𝐺𝐵(𝑣,𝑠) from (5.33) and 𝑍_𝐺𝐺(𝑣,𝑠)=𝐶_𝐺𝐺(𝑣,𝑠) from (5.34), completing the proof for (5.27).

5.7. Derivation of MAP

We recall that an MAP is constructed from an absorbing Markov-chain 𝐽^∗(𝑡) in continuous time on 𝒥=𝐺∪𝐵, with 𝐵 being a set of absorbing states, governed by 𝜈^∗ defined in (2.32). A replacement Markov-chain 𝐽(𝑡) is then generated from 𝐽^∗(𝑡), where 𝐽(𝑡) coincides with 𝐽^∗(𝑡) within 𝐺 starting from a state in 𝐺. As soon as 𝐽^∗(𝑡) reaches state 𝑟∈B, it is instantaneously replaced at state 𝑗∈𝐺 with probability 𝑝̃_𝑟𝑗 and the process continues.

In order to deduce an MAP from the unified multivariate counting process [𝑀(𝑡),𝑁(𝑡)] as a special case, we start from (4.20) in Theorem 4.1, where the underlying semi-Markov process is now the replacement Marcov chain 𝐽(𝑡) discussed above. This replacement Marcov chain is defined on 𝐺 governed by 𝜈=𝐶+𝐷 with 𝐶=[𝑐_𝑖𝑗]=[𝜈_𝑖𝑗]_{𝑖,𝑗∈𝐺} and 𝐷=[𝑑_𝑖𝑗]=[∑_𝑟∈𝐵𝜈_𝑖𝑟𝑝̃_𝑟𝑗]_{𝑖,𝑗∈𝐺} as in (2.33). We note that 𝛽_𝑖𝑗(1,𝑠)=𝛼_𝑖𝑗(𝑠) from (3.7), (4.6), (4.9) and (4.10), and 𝛽^∗𝑖(1,𝑠)={1−𝛼_𝑖(𝑠)}/𝑠 from (3.7), (4.7), (4.11) and (4.12). Hence, it follows that 𝜑1,𝑣,0+,𝑠=𝐼−𝑣_𝑖𝑗⋅𝛼_𝑖𝑗(𝑠)⁻¹𝛿^{𝑖=𝑗}⋅1−𝛼^𝑖(𝑠)𝑠.(5.44) Let 𝜈_𝐷=𝐶_𝐷+𝐷_𝐷 as in (2.35). Since 𝐽(𝑡) is a Markov chain, the dwell time in state 𝑖 is independent of the next state and is exponentially distributed with parameter 𝜈_𝑖=𝑐_𝑖+𝑑_𝑖. Consequently, one has 𝛼_𝑖(𝑠)=𝜈^𝑖𝑠+𝜈_𝑖;𝛼_𝑖𝑗(𝑠)=𝑐^𝑖𝑗+𝑑_𝑖𝑗𝜈_𝑖𝛼_𝑖(𝑠)=𝑐^𝑖𝑗+𝑑_𝑖𝑗𝑠+𝜈_𝑖.(5.45) Substituting (5.45) into (5.44) and noting [𝛿_{𝑖=𝑗}⋅((1−𝛼_𝑖(𝑠))/𝑠)]⁻¹=𝑠𝐼+𝜈_𝐷, it follows that 𝜑1,𝑣,0+,𝑠=𝑠𝐼−𝑄+𝐶+𝐷−𝑣_𝑖𝑗(𝑐_𝑖𝑗+𝑑_𝑖𝑗)⁻¹,(5.46) where 𝑄 in (2.34) is employed.

As it stands, the Laplace transform generating function of (5.46) describes the matrix counting process 𝑁(𝑡)=[𝑁_𝑖𝑗(𝑡)] where 𝑣_𝑖𝑗 corresponds to 𝑁_𝑖𝑗(𝑡). For 𝑁_MAP(𝑡) of Section 2.7, it is only necessary to count the number of replacements in (0,𝑡]. Given that state 𝑗∈𝐺 is visited from the current state 𝑖∈𝐺, this move is direct within 𝐺 with probability 𝑐_𝑖𝑗/(𝑐_𝑖𝑗+𝑑_𝑖𝑗), and such a move involves replacement with probability 𝑑_𝑖𝑗/(𝑐_𝑖𝑗+𝑑_𝑖𝑗). Accordingly, we set 𝑣_𝑖𝑗=𝑐^𝑖𝑗𝑐_𝑖𝑗+𝑑_𝑖𝑗+𝑑^𝑖𝑗𝑐_𝑖𝑗+𝑑_𝑖𝑗𝑣.(5.47) Substitution of (5.47) into (5.46) then leads to 𝜑1,𝑣,0+,𝑠=𝑠𝐼−𝑄+(1−𝑣)𝐷⁻¹,(5.48) which coincides with 𝜑(𝑣,𝑠) of (2.44), as expected.

5.8. Derivation of ACPGRP

The age-dependent counting process of Sumita and Shanthikumar [12] has been extended in this paper where the underlying renewal process is replaced by a semi-Markov process with state dependent nonhomogeneous hazard functions. Accordingly, the original model can be treated as a special case of the unified multivariate counting process proposed in this paper by setting 𝒥={0}, 𝑀(t)=𝑀₀(𝑡), 𝑁(𝑡)=𝑁₀₀(𝑡), 𝑋(𝑡)=𝑋₀(𝑡), 𝑢₀=𝑢, 𝑣₀₀=𝑣. With this specification, from Theorem 4.1, one sees that 𝜑₀₀(𝑢,𝑣,𝑤,𝑠)=𝛽^∗(𝑢,𝑤+𝑠)1−𝑣𝛽(𝑢,𝑠).(5.49) It then follows that 𝜋(𝑢,𝑠)=𝜑₀₀(𝑢,1,0,𝑠)=𝛽^∗(𝑢,𝑠)1−𝛽(u,𝑠),(5.50) which coincides with Equation (2.53).

6. Asymptotic Analysis

Let 𝒜,ℳ and 𝒩 be arbitrary subsets of the state space 𝒥 of the underlying semi-Markov process, and define 𝑀_𝒜(𝑡)=_𝑖∈𝒜𝑀_𝑖(𝑡);𝑁_ℳ𝒩(𝑡)=_{𝑖∈ℳ𝑗}_∈𝒩𝑁_𝑖𝑗(𝑡),(6.1) where 𝑀_𝒜(𝑡) describes the total number of items arrived in [0,𝑡] according to the nonhomogeneous Poisson processes within 𝒜 and 𝑁_ℳ𝒩(𝑡) denotes the number of transitions from any state in ℳ to any state in 𝒩 occurred in [0,𝑡]. Appropriate choice of 𝒜,ℳ and 𝒩 would then enable one to analyze processes of interest in a variety of applications. In the variable bit rate coding scheme for video transmission explained in Section 1, for example, one may be interested in the packet arrival stream for a specified mode of the encoder represented by 𝑀_𝒜(𝑡). In a reliability model, the underlying semi-Markov process may describe the state of a production machine. Minimal repair would take place whenever the system state is in 𝒜 at the cost of 𝑐, while complete overhaul would be forced at the cost of 𝑑 if the machine state enters 𝒩⊂𝒥. The total maintenance cost 𝑆(𝑡) is then given by 𝑆(𝑡)=𝑐𝑀_𝒜(𝑡)+𝑑𝑁_𝒥𝒩(𝑡). A simplified version of this cost structure has been analyzed by Sumita and Shanthikumar [12] where the underlying semi-Markov process is merely a renewal process with 𝒥=𝒜=𝒩={1}. The purpose of this section is to study a more general cost structure specified by 𝑆(𝑡)=𝑐𝑀_𝒜(𝑡)+𝑑𝑁_ℳ𝒩(𝑡),(6.2) with focus on the Laplace transform generating function and the related moment asymptotic behaviors of 𝑀_𝒜(𝑡),𝑁_ℳ𝒩(𝑡) and 𝑆(𝑡) based on Theorem 4.1.

For notational simplicity, we introduce the following vectors and matrices. Let 𝒜,ℳ and 𝒩⊂𝒥 with their compliments defined respectively by 𝒜^𝐶=𝒥⧵𝒜,ℳ^𝐶=𝒥⧵ℳ and 𝒩^𝐶=𝒥⧵𝒩. The cardinality of a set 𝐴 is denoted by ∣𝐴∣. 1𝑢=⎡⎢⎢⎣⋮11⎤⎥⎥⎦∈ℝ^𝐽+1;1=⎡⎢⎢⎣⋮11⋱⋯⋯⋮11⎤⎥⎥⎦∈ℝ^{(𝐽+1)×(𝐽+1)},(6.3)𝑣(𝒜)=𝑢_𝑖∈ℝ^𝐽+1with𝑢_𝑖=𝑢1,,iefls𝑖e∈;𝒜,(ℳ,𝒩)=𝑣_𝑖𝑗∈ℝ^{(𝐽+1)×(𝐽+1)}with𝑣_𝑖𝑗=𝑣1,,iefls𝑖e∈.ℳ,𝑗∈𝒩,(6.4) Submatrices of 𝐴∈ℝ^{(𝐽+1)×(𝐽+1)} are denoted by 𝐴𝐴_{𝒜••𝒜}=𝐴_𝑖𝑗_{𝑖∈𝒜,𝑗∈𝒥}∈ℝ^{∣𝒜∣×(𝐽+1)};𝐴_ℳ𝒩=𝐴_𝑖𝑗_{𝑖∈𝒥,𝑗∈𝒜}∈ℝ^{(𝐽+1)×∣𝒜∣};=𝐴_𝑖𝑗_{𝑖∈ℳ,𝑗∈𝒩}∈ℝ^{∣ℳ∣×∣𝒩∣},(6.5) so that one has

𝐴=𝐴^ℳ𝒩𝐴_ℳ𝒩𝐶𝐴_ℳ𝐶𝒩𝐴_{ℳ𝐶𝒩𝐶},(6.6) with understanding that the states are arranged appropriately.

Let 𝐴_𝑘=∫^∞0𝑥^𝑘𝑎(𝑥)𝑑𝑥, 𝑘=0,1,2,…. Throughout the paper, we assume that ‖𝐴_𝑘‖<∞ for 0≤𝑘≤2. In particular, one has 𝐴₀=𝐴(∞) which is stochastic. The Taylor expansion of the Laplace transform 𝛼(𝑠) is then given by 𝛼(𝑠)=𝐴₀−𝑠𝐴₁+𝑠²2𝐴²+𝑜𝑠².(6.7) Let 𝑒^⊤ be the normalized left eigenvector of 𝐴₀ associated with eigenvalue 1 so that 𝑒^⊤𝐴₀=𝑒^⊤ and 𝑒^⊤1=1.

We recall from Theorem 4.1 that 𝜑𝑢,𝑣,𝑤,𝑠=𝐼−𝛽̃𝑢,𝑣,𝑠⁻¹𝛽^∗𝐷𝑢,𝑤+𝑠,(6.8) where 𝛽̃𝑢,𝑣,𝑠=𝑣_𝑖𝑗⋅𝛽_𝑖𝑗𝑢_𝑖,𝑠.(6.9) From (3.6), (3.7), (4.6), (4.9) and (4.10), one sees that 𝛽_𝑖𝑗(𝑢,𝑠)=^∞0𝑒^−𝑠𝑡𝑒^{−𝐿𝑖(𝑡)(1−𝑢)}𝑎_𝑖𝑗(𝑡)𝑑𝑡.(6.10) Similarly, it follows from (3.6), (3.7), (4.7), (4.11) and (4.12) that 𝛽^∗𝑖(𝑢,𝑠)=^∞0𝑒^−𝑠𝑡𝑒^{−𝐿𝑖(𝑡)(1−𝑢)}𝐴_𝑖(𝑡)𝑑𝑡.(6.11) The 𝑟th order partial derivatives of 𝛽_𝑖𝑗(𝑢,𝑠) and 𝛽^∗𝑖(𝑢,𝑠) with respect to 𝑢 at 𝑢=1 are then given by 𝜁_{𝑟∶𝑖𝑗}(𝑠)=^def𝜕𝜕𝑢^𝑟𝜁𝛽_𝑖𝑗(𝑢,𝑠)|||_𝑢=1=^∞0𝑒^−𝑠𝑡𝐿^𝑟𝑖(𝑡)𝑎_𝑖𝑗(𝑡)𝑑𝑡,𝑟=1,2;_{∗𝑟∶𝑖}(𝑠)=^def𝜕𝜕𝑢^𝑟𝛽^∗𝑖(𝑢,𝑠)|||_𝑢=1=^∞0𝑒^−𝑠𝑡𝐿^𝑟𝑖(𝑡)𝐴_𝑖(𝑡)𝑑𝑡,𝑟=1,2.(6.12) In matrix form, (6.12) can be written as 𝜁𝜁_𝑟(𝑠)=^∞0𝑒^−𝑠𝑡𝐿^𝑟𝐷(𝑡)𝑎(𝑡)𝑑𝑡,𝑟=1,2;_{∗𝑟∶𝐷}(𝑠)=^∞0𝑒^−𝑠𝑡𝐿^𝑟𝐷(𝑡)𝐴_𝐷(𝑡)𝑑𝑡,𝑟=1,2.(6.13)

Let Φ_𝑟∶𝑘=∫^∞0𝑡^𝑘𝐿^𝑟𝐷(𝑡)𝑎(𝑡)𝑑𝑡 and Φ^{∗𝑟∶𝐷∶𝑘}=∫^∞0𝑡^𝑘𝐿^𝑟𝐷(𝑡)𝐴_𝐷(𝑡)𝑑𝑡, 𝑟=1,2. The Taylor expansion of 𝜁_𝑟(𝑠) and 𝜁^{∗𝑟∶𝐷}(𝑠) is then given by 𝜁_𝑟(𝑠)=Φ_𝑟∶0−𝑠Φ_𝑟∶1+𝑠²2Φ^𝑟∶2𝜁+𝑜𝑠²,𝑟=1,2;_{∗𝑟∶𝐷}(𝑠)=Φ^{∗𝑟∶𝐷∶0}−𝑠Φ^{∗𝑟∶𝐷∶1}+𝑠²2Φ^{∗𝑟∶𝐷∶2}+𝑜𝑠²,𝑟=1,2.(6.14)

In order to prove the main results of this section, the following theorem of Keilson [23] plays a key role.

Theorem 6.1 (Keilson [23]). As 𝑠→0+, one has 𝐼−𝛼(𝑠)⁻¹=1𝑠𝐻¹+𝐻₀+𝑜(1),(6.15) where 𝐻₁=1𝑚1𝐻₀𝑒^⊤,𝑚=𝑒^⊤𝐴₁1,=𝐻₁−𝐴₁+12𝐴²𝑍₀𝐻₁+𝑍₀−𝐻₁𝐴₁𝑍₀𝐴₀−𝐴₁𝐻₁+𝐼,=𝐼−𝐴₀+1⋅𝑒^⊤⁻¹.(6.16)

We are now in a position to prove the first key theorem of this section.

Theorem. Let 𝑝^⊤(0) be an initial probability vector of the underlying semi-Markov process. As 𝑡→∞, one has EE𝑁𝑀_𝒜(𝑡)=𝑝^⊤(0)𝑡𝑃₁+𝑃₀1+𝑜(1),_ℳ𝒩(𝑡)=𝑝^⊤(0)𝑡𝑄₁+𝑄₀1+𝑜(1),(6.17) where 𝑃₁=𝐻_1∶•𝒜Φ_{1∶0∶𝒜•},𝑃𝑄₀₁=𝐻_0∶•𝒜Φ_{1∶0∶𝒜•}−𝐻_1∶•𝒜Φ_{1∶1∶𝒜•}+𝐻_1∶•𝒜Φ^{∗1∶𝐷∶0∶𝒜•},𝑄₀=𝐻_1∶•ℳ𝐴_0∶ℳ𝒩,0_ℳ𝒩𝐶,=𝐻_0∶•ℳ𝐴_0∶ℳ𝒩,0_ℳ𝒩𝐶−𝐻_1∶•ℳ𝐴_1∶ℳ𝒩,0_ℳ𝒩𝐶.(6.18)

Proof. We first note from (6.8) together with (6.4) that ℒℒEE𝑣𝑢^{𝑀𝒜(𝑡)}_{𝑁ℳ𝒩(𝑡)}=𝑝^⊤(0)𝜑𝑢(𝒜),1,0,𝑠1;=𝑝^⊤(0)𝜑1,𝑣(ℳ,𝒩),0,𝑠1.(6.19) By taking the partial derivatives of (6.19) with respect to 𝑢 at 𝑢=1 and 𝑣 at 𝑣=1 respectively, one has ℒE𝑀_𝒜(𝑡)=𝑝^⊤(0)𝐼−𝛼(𝑠)⁻¹1𝑠𝜁^1∶𝒜•0(𝑠)_𝒜𝐶•+𝜁^{∗1,𝐷∶𝒜•}0(𝑠)_𝒜𝐶•1ℒE𝑁_ℳ𝒩;(𝑡)=1𝑠𝑝^⊤(0)𝐼−𝛼(𝑠)⁻¹𝛼^ℳ𝒩(𝑠)0_ℳ𝒩𝐶0_ℳ𝐶𝒩0_{ℳ𝐶𝒩𝐶}1.(6.20) Theorem 6.1 of Keilson [23] combined with (6.7), (6.13) then yields the Laplace transform expansions of (6.20), and the theorem follows by taking the inversion of the Laplace transform expansions.The next theorem can be shown in a similar manner by differentiating (6.19) twice with respect to 𝑢 at 𝑢=1 and 𝑣 at 𝑣=1 respectively, and proof is omitted.Theorem. As 𝑡→∞, one has EE𝑀𝑁_𝒜(𝑡)𝑀_𝒜(𝑡)−1=𝑝^⊤(0)𝑡²𝑃²¹+𝑡2𝑃₁𝑃₀+2𝑃₀𝑃₁+𝑃₂1+𝑜(𝑡);_ℳ𝒩(𝑡)𝑁_ℳ𝒩(𝑡)−1=𝑝^⊤(0)𝑡²𝑄²¹+2𝑡𝑄₁𝑄₀+𝑄₀𝑄₁1+𝑜(𝑡),(6.21) where 𝑃₀=𝐻_0∶•𝒜Φ_{1∶0∶𝒜•}−𝐻_1∶•𝒜Φ_{1∶1∶𝒜•}, 𝑃₂=𝐻_1∶•𝒜Φ_{2∶0∶𝒜•} and other matrices are as defined in Theorem 6.2.Theorems 6.2 and 6.3 then lead to the following theorem providing the asymptotic expansions of Var[𝑀_𝒜(𝑡)] and Var[𝑁_ℳ𝒩(𝑡)].Theorem. As 𝑡→∞, one has VVaarr𝑁𝑀_𝒜ℳ𝒩(𝑡)=𝑡𝑝^⊤(0)𝑈₀1+𝑜(𝑡),(𝑡)=𝑡𝑝^⊤(0)𝑉₀1+𝑜(𝑡),(6.22) where 𝑈₀=2𝑃𝑉₁₀𝑃₀+𝑃₁−𝑃₁1⋅𝑝^⊤(0)𝑃₀+2𝑃₀𝑃₁−𝑃₀𝑃₁+𝑃₂;=2𝑄₁𝑄₀+𝑄₁−𝑄₁1⋅𝑝^⊤(0)𝑄₀+𝑄₀𝑄₁.(6.23)

Proof. It can be readily seen that Var[𝑋]=E[𝑋]²−E[𝑋]²=E[𝑋(𝑋−1)]+E[𝑋]−E[𝑋]².(6.24) Substituting the results of Theorems 6.2 and 6.3 into this equation, one sees that VVaarr𝑁𝑀_𝒜ℳ𝒩(𝑡)=𝑝^⊤(0)𝑡²𝑈₁+𝑡𝑈₀1+𝑜(𝑡),(𝑡)=𝑝^⊤(0)𝑡²𝑉₁+𝑡𝑉₀1+𝑜(𝑡),(6.25) where 𝑈₀𝑈₁=𝑃₁𝑃₁−1⋅𝑝^⊤(0)𝑃₁,=2𝑃₁𝑃₀+𝑃₁𝑉−₁P₁1⋅𝑝^⊤(0)𝑃₀+2𝑃₀𝑃₁−𝑃₀1⋅𝑝^⊤(0)𝑃₁+𝑃₂;=𝑄₁𝑄₁−1⋅𝑝^⊤(0)𝑄₁𝑉₀,=2𝑄₁𝑄₀+𝑄₁−𝑄₁1⋅𝑝^⊤(0)𝑄₀+2𝑄₀𝑄₁−𝑄₀1⋅𝑝^⊤(0)𝑄₁.(6.26) From Theorem 6.1, one has 𝐻₁=(1/𝑚)1𝑒^⊤ which is of rank one having identical rows. As can be seen from Theorem 6.2, both 𝑃₁ and 𝑄₁ satisfy 1⋅𝑝^⊤(0)𝑃₁=𝑃₁,1⋅𝑝^⊤(0)𝑄₁=𝑄₁,(6.27) so that 𝑈₁=𝑉₁=0 and the theorem follows.The asymptotic behavior of E[𝑆(𝑡)] can be easily found from (6.2) and Theorem 6.2. The asymptotic expansion of Var[𝑆(𝑡)], however, requires a little precaution because it involves the joint expectation of 𝑀_𝒜(𝑡) and 𝑁_ℳ𝒩(𝑡). More specifically, one has Var[𝑆(𝑡)]=E𝑆(𝑡)²−E[𝑆(𝑡)]²=E𝑐𝑀_𝒜(𝑡)+𝑑𝑁_ℳ𝒩(𝑡)²−E𝑐𝑀^𝒜(𝑡)+𝑑𝑁_ℳ𝒩(𝑡)²=−𝑐²E𝑐₂𝑀^2𝒜(𝑡)+2𝑐𝑑E𝑀_𝒜(𝑡)𝑁_ℳ𝒩(𝑡)+𝑑²E𝑁^2ℳ𝒩(𝑡)E𝑀_𝒜(𝑡)²−2𝑐𝑑E𝑀_𝒜(𝑡)E𝑁_ℳ𝒩(𝑡)−𝑑²E𝑁_ℳ𝒩(𝑡)²,(6.28) so that Var[𝑆(𝑡)]=+𝑐²2V𝑐𝑑arE𝑀𝑀_𝒜𝒜(𝑡)+𝑑²Var𝑁_ℳ𝒩(𝑡)(t)𝑁_ℳ𝒩(𝑡)−2𝑐𝑑E𝑀_𝒜(𝑡)E𝑁_ℳ𝒩(𝑡).(6.29) In order to evaluate E[𝑀_𝒜(𝑡)𝑁_ℳ𝒩(𝑡)], we note from (6.8) that ℒE𝑀_𝒜(𝑡)𝑁_ℳ𝒩(𝑡)=𝑝^⊤(0)𝜕²𝜕𝑢𝜕𝑣𝜑(𝑢(𝒜),𝑣(ℳ,𝒩),0+,𝑠)|||_𝑢=𝑣=11.(6.30) The asymptotic expansion of E[𝑀_𝒜(𝑡)𝑁_ℳ𝒩(𝑡)] can then be obtained as for the previous theorems.Theorem. As 𝑡→∞, one has E𝑀_𝒜(𝑡)𝑁_ℳ𝒩(𝑡)=𝑝^⊤(0)𝑡²2𝑇¹+𝑡𝑇₀1+𝑜(𝑡),(6.31) where 𝑇₁=𝑃₁𝑄₁𝑅+𝑄₁𝑃₁,𝑇₀=𝑃₀𝑄₁+𝑃₁𝑄₀+𝑅+𝑄₀𝑃₁+𝑄₁𝑃₀,=𝐻_{1∶•,(𝒜∩ℳ)}Φ^{1∶0∶(𝒜∩ℳ),𝒩},0_{(𝒜∩ℳ),𝒩𝐶}.(6.32)Now the key theorem of this section is given from (6.29), Theorems 6.2, 6.4 and 6.5.Theorem. As 𝑡→∞, one has E[𝑆(𝑡)]=𝑝^⊤(0)𝑡𝑐𝑃₁+𝑑𝑄₁+𝑐𝑃₀+𝑑𝑄₀1Var[𝑆(𝑡)]=𝑡𝑝_⊤+𝑜(1),(6.33)(0)𝑊₀1+𝑜(𝑡),(6.34) where 𝑊₀=𝑐²𝑈₀+𝑑²𝑉₀+2𝑐𝑑𝑇₀−2𝑐𝑑(𝑃₁1⋅𝑝^⊤(0)𝑄₀+𝑃₀1⋅𝑝^⊤(0)𝑄₁).Proof. Equation (6.33) follows trivially from Theorem 6.2. For (6.34), we note from Theorems 6.2, 6.4 and 6.5 together with (6.29) that Var[𝑆(𝑡)]=𝑝^⊤(0)𝑡²𝑊₁+𝑡𝑊₀1+𝑜(𝑡),(6.35) where 𝑊₀𝑊₁=𝑐𝑑𝑇₁−2𝑐𝑑𝑃₁1⋅𝑝^⊤(0)𝑄₁,=𝑐²𝑈₀+𝑑²𝑉₀+2𝑐𝑑𝑇₀−2𝑐𝑑𝑃₁1⋅𝑝^⊤(0)𝑄₀+𝑃₀1⋅𝑝^⊤(0)𝑄₁.(6.36) From (6.27), the first term on the right-hand side of Equation (6.35) can be rewritten as 𝑐𝑑𝑝^⊤(0)𝑇₁1−2𝑐𝑑𝑝^⊤(0)𝑃₁1⋅𝑝^⊤(0)𝑄₁1=𝑐𝑑𝑝^⊤(0)𝑃₁𝑄₁+𝑄₁𝑃₁1−2𝑐𝑑𝑝^⊤(0)𝑃₁1𝑝^⊤(0)𝑄₁1=𝑐𝑑𝑝_⊤(0)𝑃₁𝑄₁1+𝑝^⊤(0)𝑄₁𝑃₁1−2𝑐𝑑𝑝^⊤(0)𝑃₁1𝑝^⊤(0)𝑄₁1=𝑐𝑑𝑝_⊤(0)𝑃₁1𝑝^⊤(0)𝑄₁1+𝑝^⊤(0)𝑄₁1𝑝^⊤−2𝑐𝑑𝑝_⊤(0)𝑃₁1(0)𝑃₁1𝑝^⊤(0)𝑄₁1=0,(6.37) completing the proof.

7. Dynamic Analysis of a Manufacturing System for Determining Optimal Maintenance Policy

As an application of the unified multivariate counting process, in this section, we consider a manufacturing system with a certain maintenance policy, where the system starts anew at time 𝑡=0, and tends to generate product defects more often as time goes by. When the system reaches a certain state, the manufacturing system would be overhauled completely and the system returns to the fresh state. More specifically, let 𝐽(𝑡) be a semi-Markov process on 𝒥={0,1,2,…,𝐽} governed by 𝐴(𝑥), describing the system state at time 𝑡 where state 0 is the fresh state and state 𝐽 is the maintenance state. When the system is in state 𝑗, 0≤𝑗≤𝐽−1, product defects are generated according to an NHPP with intensity 𝜆_𝑗(𝑥). It is assumed that the system deteriorates monotonically and accordingly 𝜆_𝑗(𝑥) increases as a function of both 𝑥 and 𝑗. When the system reaches state 𝐽, the manufacturing operation is stopped and the system is overhauled completely. The maintenance time increases stochastically as a function of 𝐽. In other words, the further the maintenance is delayed, the longer the maintenance time would tend to be. Upon finishing the overhaul, the system is brought back to the fresh state 0. Of interest, then, is to determine the optimal maintenance policy concerning how to set 𝐽.

In order to determine the optimal maintenance policy, it is necessary to define the objective function precisely. Let 𝜓_𝑑 be the cost associated with each defect and let 𝜓_𝑚 be the cost for each of the maintenance operation. If we define two counting processes 𝑀(𝑡) and 𝑁(𝑡) as the total number of defects generated by time 𝑡 and the number of the maintenance operations occurred by time 𝑡 respectively, the total cost in [0,𝑇] can be described as 𝐶_𝐽(𝑇)=𝜓_𝑑⋅E[𝑀(𝑇)]+𝜓_𝑚⋅E[𝑁(𝑇)].(7.1) Let 𝑁 be the set of natural numbers. The optimal maintenance policy 𝐽^∗ is then determined by 𝐶_𝐽^∗(𝑇)=min_𝐽∈𝑁𝐶_𝐽(𝑇).(7.2)

In what follows, we present a numerical example by letting 𝒥={0,1,…,𝐽} for 1≤𝐽≤9. For the underlying semi-Markov process, we define the matrix Laplace transform 𝛼(𝑠) having IFR (Increasing Failure Rate) and DFR (Decreasing Failure Rate) dwell time distributions as described below, where the underlying parameters are set in such a way that the means of IFR dwell times are equal to those of DFR dwell times. By introducing matrices 𝜃,𝜃̂ and 𝑝, for which the details are given in Table 1 along with other parameter values, we define 𝛼_IFR(𝑠)=𝜃^𝑖𝑖𝑠+𝜃_𝑖𝑖⋅𝜃^𝑖𝑗𝛼_DFR𝑠+∑^𝐽𝑗=0𝜃_𝑖𝑗,(𝑠)=𝑝_𝑖⋅𝜃^𝑖𝑗𝑠+∑^𝐽𝑗=0𝜃_𝑖𝑗+1−𝑝_𝑖⋅𝜃̂^𝑖𝑗𝑠+∑^𝐽𝑗=0𝜃̂_𝑖𝑗,(7.3) where 𝑝_𝑖=𝜃^𝑖𝑖+∑^𝐽𝑗=0𝜃_𝑖𝑗/𝜃^𝑖𝑖⋅∑^𝐽𝑗=0𝜃_𝑖𝑗−1/∑^𝐽𝑗=0𝜃̂_𝑖𝑗1/∑^𝐽𝑗=0𝜃_𝑖𝑗−1/∑^𝐽𝑗=0𝜃̂_𝑖𝑗,^𝐽_𝑗=0𝜃_𝑖𝑗≠^𝐽_𝑗=0𝜃̂_𝑖𝑗.(7.4)

Table 1 

Parameters of the numerical example.

Based on Theorem 6.2, the asymptotic behaviors of the mean of 𝑀(𝑡) and 𝑁(𝑡) per unit time with maintenance policy 𝐽=1,…,9 are depicted in Figures 3 and 4. One could see that both the mean of 𝑀(𝑡) and 𝑁(𝑡) per unit time converges to a positive value as time 𝑡 increases. In order to determine the optimal maintenance policy, for 𝐽∈{1,2,…,9}, the corresponding total cost 𝐶_𝐽(𝑇) can be computed based on Theorem 6.6. Numerical results are shown in Table 2 and depicted in Figure 5. For the case of IFR, the optimal maintenance policy is at 𝐽^∗=6, while 𝐽^∗=5 for the DFR case, where the running period 𝑇 is taken to be 𝑇=1000 hours.

Table 2 

Results of the problem of determining optimal maintenance policy.

(a)

(b)

Figure 3 

Mean of M(t) per unit time.

(a)

(b)

Figure 4 

Mean of N(t) per unit time.

(a)

(b)

Figure 5 

Optimal maintenance policy: IFR versus DFR with T = 1000.

8. Concluding Remarks

In this paper, a unified multivariate counting process [𝑀(𝑡),𝑁(𝑡)] is proposed with nonhomogeneous Poisson processes lying on a finite semi-Markov process. Here the vector process 𝑀(𝑡) counts the cumulative number of such nonhomogeneous Poisson arrivals at every state and the matrix process 𝑁(𝑡) counts the cumulative number of state transitions of the semi-Markov process in [0,𝑡]. This unified multivariate counting process contains many existing counting processes as special cases. The dynamic analysis of the unified multivariate counting process is given, demonstrating the fact that the existing counting processes can be treated as special cases of the unified multivariate counting process. The asymptotic behaviors of the mean and the variance of the unified multivariate counting process are analyzed. As an application, a manufacturing system with certain maintenance policies is considered. The unified multivariate counting process enables one to determine the optimal maintenance policy minimizing the total cost. Numerical examples are given with IFR and DFR dwell times of the underlying semi-Markov process. As for the future agenda, the impact of such distributional properties on the optimal maintenance policy would be pursued theoretically. Other possible theoretical extensions include: (1) analysis of the reward process associated with the unified multivariate counting process; and (2) exploration of further applications in the areas of modern communication networks and credit risk analysis such as CDOs (collateralized debt obligations) for financial engineering.

Acknowledgment

The authors wish to thank two anonymous referees and the associate editor for providing constructive comments and valuable references, which contributed to improve the original version of the paper. This research is supported by MEXT Grand-in-Aid for Scientific Research (C) 17510114.