NOTES ON MARKOVIAN POINT
PROCESSES
Introduction
Renewal processes provide simple models of point processes which may
describe an ordered set of points (arrival instants, service completion
epochs, equipment failures, etc.) on [0, oo). Their main simplifying
feature is the independence and equidistribution of successive interrenewal
intervals. We study in some detail renewal processes for which
the intervals have a PH distribution.
In the second part of the chapter, we describe a Markovian point
process which generalizes and unifies the Markovian processes commonly
used in applied probability.
The results presented here illustrate anew the simplicity and the
algorithmic tractability offered by the underlying Markovian structure
in the method of phases.
3.1 The PH Renewal Process
We consider a renewal process with a PH distribution for the interrenewal
intervals. We call such a process a PH renewal process. We shall
primarily deal with the continuous case, stating the parallel results for
discrete time only in summary form.
As with PH distributions, we may associate a Markov process with a
PH renewal process. Assume as given a Markov process on {0,1,..., n}
61
62 Chapter 3. Markovian Point Processes
with initial distribution (0, T) and infinitesimal generator
which is associated with the PH distribution F(-); we assume that
there is no atom at 0. Start the process at time 0 and let it evolve
to absorption. Instantaneously upon absorption, restart the Markov
process by choosing a new state with the distribution r and let it
proceed to absorption again, etc.
Consider the set {0 = to < ti < t-2 < • • •} of time points at which
the Markov process is reinitialized with the distribution T. This clearly
forms a renewal process with interrenewal distribution PH(r,jT).
If the state at time 0 is chosen with a different distribution j3, then
we have a delayed renewal process where the first interval has the distribution
PH(/3, T). Unless otherwise stated, when we refer to a renewal
process in what follows, we consider the nondelayed case, i.e., the one
where there is a renewal at time 0.
With respect to renewal processes, some key quantities of interest
are the renewal function R(x), which is the expected number of renewals
in the interval [0,a;], and the renewal density r(x) = R'(x), which has
the physical interpretation that r(x}dx is the elementary probability of
having a renewal in the interval (x, x + dx).
Let us begin by observing that our construction above defines a
Markov process {J(t) : t > 0} on {!,...,n}; we ignore the instantaneous
visits to the state 0 and require that the trajectories be right
continuous. This is the "phase process" associated with the PH renewal
process. Its infinitesimal generator is given by the matrix
We see this by noting that in our construction of the PH renewal process,
there may be a transition from i to j with i ^ j in two ways: either
directly with the rate T^- or after absorption in the state 0, which occurs
at the rate ^ instantaneously followed by a restart in the state j,
which has probability TJ. The diagonal elements of D are such that
DI = 0. This Markov process plays a key role in analyzing the PH
renewal process.
3.1. The PH Renewal Process 63
Recall that a phase j is transient if there is a path from j to the
absorbing phase 0. Symmetrically, the phase j is said to be useful if
there exists a path to j from some phase i with TI > 0; if there were
no such i, then j would never be visited and it would be useless to
include it in the model. It will be much simpler, from now on, to
restrict ourselves to irreducible PH representations, i.e., to pairs (r,T)
such that all phases are transient and all are useful. Clearly, there is
no loss of generality if we only consider irreducible representations (see
Theorem 2.4.2). The following property is immediate.
Lemma 3.1.1 If the representation (T, T) is irreducible, then the matrix
D = T + t • T is irreducible.
The renewal density of a PH renewal process is expressed as follows.
Theorem 3.1.2 Consider the (possibly delayed) PH renewal process
with representation (T, T) and initial phase distribution (3. Its renewal
density is given by r(x) = (3exp(Dx)t, x > 0.
Proof. A renewal occurs in (x, x+dx) if and only if the phase process
{J ( t ) : t > 0} is in some state j at time x, and an absorption into 0
occurs in (x, x + dx). Therefore,
Setting (3 = T gives the renewal density of the (nondelayed) renewal
process which has a renewal at time 0. With this, we may determine
the renewal function (of the nondelayed renewal process) as
where r(u) = rexp(Du)t. In order to evaluate this, we shall need the
following algebraic property.
64 Chapter 3. Markovian Point Processes
Lemma 3.1.3 The matrix D — II is nonsingular, where
and
The vector TT is the stationary probability vector of the phase process,
i.e., the solution of the system TcD = 0, TT! = 1. Furthermore, we have
that 7t(D - H)-1 = -TT.
Proof. The proof that the vector TT given in (3.1) is the solution of
7TJ9 = 0, TT! = 1 is by direct verification.
Assume that there exists a nonzero vector x such that x(D— I-TT) =
0. After postmultiplication by 1, we find that xl = 0 and therefore
that xD = 0. Since D is irreducible by Lemma 3.1.1, this implies that
x is proportional to TT, which contradicts the earlier conclusion that
xl = 0; thus the first statement is proved.
Finally, TT may be written as TT = (7rl)7r — irD, or TT = 7r(II — D),
from which we have the last statement.
With this, we easily determine an expression for the renewal function.
Theorem 3.1.4 The renewal function of the PH renewal process (starting
with a renewal at time 0) with representation (r, T) is given by
where M = r(—T)"1! is the mean interrenewal time and U — 1 • TT.
Proof. Since R(x] is the expected number of renewals in [0, x] under
the assumption that the first renewal occurs at time 0, we have that
3.2. The Number of Renewals 65
If we write t as (D-U)(D-U)-lt, using the facts that DnU = Dnl-7T =
0 for n > 1 and that it(D — H)"1 = —TT, we obtain that
which concludes the proof since
by (2.13).
3.2 The Number of Renewals
Consider the PH renewal process defined by (r,T), and let N(x) denote
the number of renewals in (0, x]. Our construction clearly shows
that (-/V(x), J(x)) is Markovian. Thus, to determine the distribution of
N(x), it is convenient to consider the joint distributions
for k > 0, 1 < i,j < n. We shall denote by P(k,x) the n x n matrix
{Pij(k, x)}, and we also define the matrix generating function
The marginal distribution of N(x) is given by the scalar sequence
{(3P(k, x)l : k > 0} if (3 is the initial phase distribution.
Theorem 3.2.1 The matrices P(-,-) satisfy the following system of
differential equations:
with
Their generating function is given by
66 Chapter 3. Markovian Point Processes
Proof. The generator of the process {(N(x), J(x)) : x > 0} is given
in block form by the matrix
The matrix of transition functions is given by
since we have that
for k, k' > 0, 1 < i, j < n.
The proof of (3.2, 3.3) is immediate once we write in block form
the Kolmogorov backward differential equation P'(x) = QP(x). From
there, we find that
so that (3.4) follows.
These equations provide matrix analogues of the corresponding results
for the Poisson process (take n — 1, T = —A, and r = 1).
Remark 3.2.2 We now give an alternate, more elaborate derivation
using a renewal argument.
Clearly, we have that P(0, x) = exp(Tx) since N(x) = 0 if and only
if all the phase transitions up to time x are among the transient states.
This immediately gives the differential equation for P(0, x).
3.2. The Number of Renewals 67
For k > 1, we condition on the first renewal epoch x — u and write
that
premultiplying both sides by exp(—Tx) and differentiating with respect
to x, we obtain the differential equation for P(k,x).
Multiplying by zk and adding, we get
Now, premultiply by exp(—Tx) and differentiate; this yields (3.5).
With some experience, one can write the integral equation above
without having to go through the elaborate steps of setting up the
equations for P(k,x)\ one uses a conditioning argument on the first
arrival epoch x — it, at which one "throws in" a factor z to record the
fact that there is a renewal at that epoch. Later, when we have to write
equations for matrix generating functions, we shall often employ this
device without explaining the details.
As in the preceding chapter, uniformization offers an alternative
method to compute the probabilities P(-, •). Its advantages relative to
solving the system of differential equations by numerical integration are
similar to those found earlier when we discussed the computation of the
PH density and distribution function. With both methods, the actual
computation of the probabilities is a nontrivial task, expensive both in
memory and time requirements. Fortunately, in the models that will be
discussed, many quantities of interest may be obtained without having
to compute these quantities explicitly.
Lemma 3.2.3 Define c — max{|Tij| : 1 < i < n}, P = c~lT + I, and
p = c~lt = I- PI.
The matrices P(y, x) (y > 0, x > 0) are given by
where the matrices Km^ of order n are defined as follows:
68 Chapter 3. Markovian Point Processes
for ra, z/ > 0.
Proof. Consider the uniformized version of the phase process {J(t) :
t > 0} associated with the PH renewal process. As noted in section 2.8,
this allows us to consider the Markov process as made up of a Poisson
process and a discrete time Markov chain. The Poisson process with
rate c determines when an absorbing discrete time Markov chain makes
its transitions. The transition matrix among transient states is P,
and the absorption probabilities are given by the vector p. When the
Markov chain gets absorbed, a new transient state is instantaneously
chosen with the distribution r.
Now, given that the initial phase is i, let (Km^)ij denote the conditional
probability that the Markov chain is in the transient phase j
immediately after the rath transition and that there have been v instantaneous
visits to the absorbing state during the first ra transitions.
The equations (3.7), (3.8) are then immediate, and (3.9) is obtained by
a simple conditioning argument on the first transition.
The equation (3.6) follows immediately if we condition on the number
of Poisson events in (0, x\.
Example 3.2.4 This is a continuation of Example 2.1.2; it gives one
more illustration of our general theme, that the Markovian structure in
phase processes allows for a strong similarity to the Poisson process.
Consider the PH/G/oo queue which is the infinite server queue with
PH(T,T) renewal arrivals and general service time distribution H(-).
Let Pij(n, x) denote the probability, given that at time 0 all servers are
idle and that the arrival process is in phase i, that at time x, exactly n
servers are busy and the arrival process is in phase j. Define the matrix
generating function
As in Example 2.1.2, we condition on the first arrival; as in the proof
of Theorem 3.2.1, we keep track of the phase at epochs of arrival and
at time x. This gives
3.2. The Number of Renewals 69
In the right side of the above equation, the first term corresponds to the
case where there is no arrival in (0, x]. The second term corresponds to
the case where an arrival occurs at time x — u and leaves before time
x; at time x — u a new phase is chosen with the distribution T. In the
third case, the customer who arrives at time x — u is still present at
time x and contributes a factor z.
We proceed as in Remark 3.2.2, premultiply by exp(—Tx), differentiate,
and obtain
with the initial condition P*(2,0) = 1.
Prom (3.10), one may numerically obtain the moments of the queue
length at time x. For further details, see Ramaswami and Neuts [101].
Example 3.2.5 This is the continuation of Examples 2.4.6 and 2.5.4.
We consider here a PH renewal process with representation (T, T) and
define N as the number of renewals in the interval (0, X), where X has
an exponential distribution with rate A, independent of the renewal
process. The number of renewals has a geometric distribution, with
parameter p = r51, where 5 = A(A/ — T)"1.
To see this, we compare X and the first renewal epoch, say, V, and
obtain that P[N = 0] = P[X < Y]. If X > Y, then N = I + Nlt
where N\ is the number of renewals in the remainder of the interval
(y, X). Given the memoryless property of the exponential distribution,
NI has the same distribution as N and we find that N is geometrically
distributed with parameter p = P[N — 0]. Now,
which is easily computed.
70 Chapter 3. Markovian Point Processes
3.3 The Stationary Process
Given a renewal process with interrenewal time distribution F(.), its
stationary version is obtained by imposing an initial delay with distribution
where M is the expected interval between two renewal epochs. The
distribution G(-) is the asymptotic survival distribution. That is, if
for fixed t we define the survival random variable Vt as the length of
the interval from t to the next renewal, then one has that G(x) =
lim^ooPfl/t <4
If F(') is PH(r, T), then the distribution of the survival interval Vt
is PH(/3(£), T), where (3(t) = fiexp(Dt) is the distribution of the phase
at time t for the process {J(x) : x > 0}. Knowing (by Lemma 3.1.3)
that limt_>oo 0(t) = TT, and that TT is given by (3.1), we have proved the
following theorem.
Theorem 3.3.1 If F(-) is PH(r,T), then the distribution
where M is the expected interval between two renewal epochs, is PH(TT, T),
where TT = (rT-l1)-lrT-1.
This means that the stationary version of the PH renewal process
is the delayed renewal process obtained by choosing the initial phase
according to the stationary distribution TT.
Example 3.3.2 Consider an M/PH/1 queue, i.e., an M/G/1 queue
with PH service times. A consequence of the theorem above is that
the stationary waiting time distribution W(-) is PH. If the service time
distribution is PH(r,T), then W(-) is PH((1 - p)ir,T).
To prove this, we recall the Pollachek-Khinchine formula (Cohen
[16, Page 255]) which states that
3.4. Discrete Time PH Renewal Processes 71
where w(s) and (f)(s) are, respectively, the Laplace-Stieltjes transforms
of the stationary waiting time and of the service time distributions and
p = AM is the traffic coefficient; the expression for </>(s) is given in
(2.12).
This may also be written as
where
is the Laplace-Stieltjes transform of the asymptotic survival distribution
(3.11).
Therefore, the distribution function W(-) itself is such that
i.e., it is a geometric mixture of the convolution powers of G(-). Since
G(-) is PH by Theorem 3.3.1, and since the geometric distribution is
discrete PH, W(-) is PH by Theorem 2.6.3.
3.4 Discrete Time PH Renewal Processes
We can use discrete time Markov chains very much like we used continuous
time Markov processes in section 3.1. This allows us to define discrete
time PH renewal processes with interrenewal intervals distributed
as PHd(r,T).
We assume that rl = 1, so that renewals occur singly. If rl <
1, then renewals occur in batches of random size with a geometric
distribution. This is a special case of the process which we define in
the next section.
The generator of the phase process is T* = T + t • r with t =
1 — Tl. We may, as before, assume without loss of generality that it
is irreducible. We shall also assume that it is aperiodic. Under these
conditions, the following results are easy to prove; we omit the details.
72 Chapter 3. Markovian Point Processes
The stationary probability vector of T* is TT = (l/E[X])r(I
where E[X] — r(I — T)"1! is the expected interrenewal time. The
stationary version of the renewal process is obtained by choosing the
initial phase with the vector TT.
The renewal density {r(k} : k > 0}, where r(k) is the probability of
a renewal at epoch k, is given by r(k) = r(T -+-1 - r}k~lt for k > 1.
The distribution of the remaining lifetime at k is PH(/3fc, T), where
/3Q is the initial phase distribution and (3k = /30(T + t- r)k. As k tends
to infinity, this distribution converges to PHrf(7r, T).
The density of the elapsed lifetime Yk at k is given by
for j — k and
for 0 < j < k — 1; it converges to PHd(7r, T) as k tends to infinity.
3.5 A General Markovian Point Process
Let us consider anew the process {(N(x), J(x)) : x > 0} of section 3.2,
where N(x) is the number of renewals in (0, x] and J(x) is the phase
at time x. Its infinitesimal generator is
and is very similar to the generator of the familiar Poisson process:
Thus it appears that the PH renewal process may be seen as a generalization
of the pure birth process. There, the birth process N(x) is
3.5. A General Markovian Point Process 73
modulated by the phase process J(x) whose state dictates the instantaneous
birth rates.
If we partition the state space {(n, j) : n > 0,1 < j < m} into
subsets
called levels , we see that the diagonal blocks in (3.12) are formed by the
rates at which phase changes occur while the level remains the same.
The blocks above the diagonal are formed by the rates at which a birth
occurs: the process moves from one level to the next level above, and a
change of phase may also occur. In a PH renewal process, these phase
changes occur in a very restricted way since the new phase is chosen at
each birth with the same distribution T.
We may generalize this birth process in two ways: we may allow for
multiple births so that a transition is possible in one step from each level
l(k) to some level t(k + v] with v > 1; second, we may allow phase
changes at birth epochs to occur without any restriction. We define
in this way a general Markovian counting process with infinitesimal
generator Q given below:
where DO, DI, D-z,... are matrices of order m with D^ > 0 for all k > 1
Do(i,j) > 0 for 1 < i ^ j < m, DQ(I,I) < 0 for 1 < i < m, and
£fc>0£>fcl = 0.
The only restriction which we impose is that the transition rate from
(k, i) to (k + f, j) must be independent of k: the birth rate and the size
of the increase may depend on the phase, but not on the current level.
This family of counting processes has received several names in the
literature: versatile Markovian processes in Neuts [78] and Neuts processes
in Ramaswami [92]; at present, the name MAP (Markovian arrival
process) introduced in Lucantoni, Meier-Hellstern, and Neuts [61]
is frequently used, as is the name BMAP (batch Markovian arrival pro74
Chapter 3. Markovian Point Processes
cess). We now give a few examples to illustrate the variety of models
subsumed under it as special cases.
Example 3.5.1 Batch PH Renewal Processes. Consider a PH
renewal process with representation (T, T). Assume that at the nth
renewal epoch, a batch of size Xn arrives, where the Xn's are i.i.d. random
variables with density {p^ : k > 1}. This process is a Markovian
point process with generator of the form (3.13), where D0 = T and
Dk=pkt-r(k>l).
Alternatively, we may assume that the batch size depends both on
the phase i from which absorption occurs and the phase .;' in which
the process restarts with the distribution {pk(i,j} : k > 1}. Then
Dk = Pk ° (t - T), where we denote by A o B the matrix C such that
n. —. A. D . \^ij — siijj-'ij-
Example 3.5.2 Markov Modulated Poisson Process. This is a
Poisson process in which the instantaneous arrival rate depends on
an auxiliary Markov process which serves as a random environment.
Suppose that the auxiliary process has m states and generator G and
that when the auxiliary state is i, arrivals occur at rate A^. The resulting
process is called a Markov modulated Poisson process. It belongs to
the family defined here with D0 = G - A(A), A = A(A), and Dk = 0
for k > 2, where A (A) is the diagonal matrix with the Aj values in the
diagonal.
Markov modulated Poisson processes have been used extensively,
for example, in modeling arrivals of packet streams to communication
systems in which the rate of arrivals depends on the number of sources
(terminals, video codecs, etc.) that are active. The interrupted Poisson
process is a particular case with two phases and A2 = 0; in fact, it is a
PH renewal process with hyperexponential interrenewal intervals.6
We may add the assumption that arrivals occur in batches with
environment-dependent sizes: if the auxiliary variable is i, then the
batch is of size k with probability Pi (k) (k > 1). Then-E^ = A(A)A(p(&))
for k > 1.
6This is one more indication of the fact that PH processes do not have a unique
representation.
3.6. Analysis of the General Process 75
Example 3.5.3 PH Semi-Markov Processes. Consider a semi-
Markov process {(Xn, Jn) : n > 0} with PH interval distributions:
the conditional distribution of Xn given that Jn — i is PH(Tj,Ti) and
P[Jn+i = j\Jn = i] = Pij. This is also a Markovian point process; we
leave it to the reader to write the matrices D^, k > 0.
Example 3.5.4 Thinning and Superpositions. It is elementary to
verify that the superposition of two Markovian point processes is again
a Markovian point process (use the same construction as in Theorem
2.6.4 and consider a two-dimensional phase process composed of the
individual phases). Similarly, if the points of increase are discarded with
probability p independently of each other, then the resulting thinned
process remains a Markovian point process.
Example 3.5.5 Departure Processes. This is a continuation of Example
2.3.3. Consider an M/M/c/c+K queue (Poisson arrivals, exponential
services, c servers, and a buffer of size K). The overflow process
is a Markovian point process, of which the matrix DQ is given below in
the special case where c = 3 and K = 2:
the only nonzero element of D\ is A in the lower right corner.
The departure process is also a Markovian point process, as one may
easily verify. If we assume that waiting customers become impatient
and renege independently of each other and if their waiting time exceeds
an exponentially distributed interval, then the process of reneging
customers is also a Markovian point process.
3.6 Analysis of the General Process
We show in this section that the analysis of the general Markovian
point process is very similar to that of the PH renewal process; naturally,
multiple births add an element of complexity. The basic random
76 Chapter 3. Markovian Point Processes
variables of interest are the phase J(x) at time x and -/V(z), the number
of steps up (or arrivals, or events, etc.) during the interval (0, x].
Equivalently, N(x) is the level reached at time x, assuming that the
process starts in the level 0 at time 0.
We assumed earlier that the intervals between renewals had a
nondefective distribution. As we saw with the Theorems 2.4.2 and
2.4.3, this is equivalent to the assumption that the matrix T on the
diagonal of (3.12) is nonsingular. Here, we shall require that N(x)
increases without bounds as x tends to infinity. In other words, the
process never gets stuck in a level; for every pair n' > n, starting from
any state (n, z), it will move to some state (n',j) in a finite time a.s.
This holds if and only if the matrix D0 on the diagonal of (3.13) is
nonsingular.
Second, we need to assume that the expected size of the jumps to
higher levels is finite so that the expected value of N(x) remains finite
over finite intervals. This holds if and only if the vector d = Z)jt>i kDkl
is finite. This automatically holds for the PH renewal process, naturally.
The phase process {J(x) : x > 0} is a Markov process with infinitesimal
generator D — Z)fc>o Dk. Our third assumption is that the matrix
D is irreducible. We saw in section 3.1 that, for a PH renewal process,
this is equivalent to requiring that the phases are all useful. Without
this assumption, the discussion of the limiting behavior of our process
becomes more involved.
Let us denote by A(x) (x > 0) the expected value of N(x) and
by X(x) its derivative. The instantaneous rate X(x) has the physical
interpretation that \(x)dx is the expected increase in the level during
(x, x + dx). The functions A(-) and A(-) are a generalization of the
renewal function R(-) and the renewal density r(-) studied in section 3.2.
Throughout this section, we shall assume that the initial state is
(W(0) = 0, J(0) = j) with probability $ with 01 = 1.
Theorem 3.6.1 The instantaneous rate of increase is given by
its asymptotic value is
3.6. Analysis of the General Process 77
where TT is the stationary probability vector of D.
The expected increase A(x) = E[AT(a;)] over the interval (0, x] is
given by
Proof. The proof mostly repeats the arguments of Theorems 3.1.2
and 3.1.4. To prove (3.15), we recall that (under our general assumptions)
the limit of P[J(x) = j|«/(0) = i] as x —> oo is TTJ independently
of i so that f)exp(Dx) converges to TT.
We leave it to the reader to develop computational procedures for
evaluating the intensity function. This may be done by solving a suitable
system of linear differential equations. Indeed, we see from (3.14)
that A(x) = 0y(x), where the vector y(x) is a solution of the system
y'(x) = Dy(x), y(0) = d.
Alternately, one may proceed by uniformization: (3.14) may also be
written as
where c > max \Da\ and D* = / + 1/cD.
The stationary version of the Markovian point process is obtained
by choosing 0 = TT. We then have that the marginal distribution of
the phase is TT exp(Dx) = TT at any time x and that the instantaneous
increase rate is constant: A(x) = A*; furthermore, the expected increase
over the interval (0, x] is linear: A(x) = A*x. Also, the increase probabilities
P[N(t + x) — N(x) = k] are indeed independent of x since we
have that
To conclude this chapter, we shall consider the time-dependent conditional
distribution of the Markovian point process. Define
78 Chapter 3. Markovian Point Processes
for k > 0, 1 < i,j < m; denote by P(k,x) the matrix {Pij(k, x)} and
by P*(z, x) the generating function
The following theorem is a direct generalization of Theorem 3.2.1.
Theorem 3.6.2 The matrices P(-, •) satisfy the following system of
differential equations:
fork>0 with P(0,0) = / and P(k, 0) = 0 for k > 1. The generatin
function F*(-, •) satisfies the integral equation
and is given by
where
Proof. The equation (3.16) is just a reformulation of Kolmogorov's
backward equations.
The equation (3.17) is obtained by conditioning on the first epoch
when there is a jump in level. The first term corresponds to the case
where the process remains in the level 0 during the whole interval [0, x]
The summand in the second term corresponds to the case where the
process jumps from level 0 to level k at time x — u, the epoch of the
first jump.
Premultiplying both sides of (3.17) by exp(—D0x) and differentiating
with respect to x gives
which, together with D(z, 0) = /, gives (3.18).
3.6. Analysis of the General Process 79
Considering the stationary version of the Markovian point process,
the same approach may be used to determine the joint generating function
of N(i) and N(x+i) (along with the necessary phases), from which
one obtains the auto-covariance function of the process. We refer the
reader to Neuts [82, Chapter X].
Needless to say, the material in this section may be adapted to the
discrete time case. The resulting processes have been called DMAPs
(discrete time Markovian arrival processes) by some authors.
