DESCRIPTORS FOR POINT PROCESSES BASED ON RUNS : THE MARKOVIAN ARRIVAL PROCESS

This paper is part of a broader investigation of properties of a point process that can be identified by imagining that the process is involved in a competition for the generation of runs of events. The general purpose of that methodology is to quantify the prevalence of gaps and bursts in realizations of the process. The Markovian arrival process (MAP) is highly versatile in qualitative behavior and its analysis is numerically tractable by matrix-analytic methods. It can therefore serve well as a benchmark process in that investigation. In this paper, we consider the MAP and a regular grid competing for runs of lengths at least r1 and r2, respectively. A run of length r in one of the processes is defined as a string of r successive events occurring without an intervening event in the other process.

that competition is a global measure of their relative burstiness.For several values of r 1 and r2, that fraction and several related quantities are easily estimated from actual or simulated traffic data.This is a clear advantage of such descriptors of burstiness.However, to explore the mathe- matical properties of these descriptors for familiar point processes and thereby to inspire confid- ence in their use, we examine the cases for which computable, tractable expressions for the descrip- tors can be derived.
The interest of the proposed approach can best be demonstrated by simulated examples and by using computer graphical displays.A verbal description of these examples is unavoidably lengthy and this is not the appropriate place for it.Concrete applications are to be the subject of future papers, but a brief description of the competition is given in Section 5.This article deals with the case where P1 is a Markovian arrival process (MAP), a broad generalization of the Poisson process that can be thoroughly analyzed by matrix-analytic techniques.The second process, which we think of as the reference process, is just the stationary renewal process with events that are a distance d apart.By a simple rescaling of time, we may obviously choose d = 1.We call P2, the stationary unit grid.That reference process is useful in the exploration of a single MAP, where the competition has a natural interpretation.For expositions of the essential properties of the MAP, see Lucantoni [2] or Neuts [3].Some issues in the quantification of the burstiness of point processes are discussed in Neuts [4].
The analysis of the case where P2 is the stationary unit grid is not particularly easy.As is well-known from other contexts, for continuous-time models, the equidistant events preclude a simple Markovian analysis.There are interesting special results for the MAP in the role of PI" These are useful to calibrate the run method of comparison when applied to other pairs of point processes.
As expected, there are major analytic simplifications when the reference process is Poisson.As always, the limitation in practical versatility of Poisson assumptions is severe.That makes the MAP such an appealing benchmark process.For the modest cost of performing calculations with matrix functions, we preserve much of the tractability of the Poisson process while retaining great modeling versatility.
In subsequent papers, we shall study the MAP with two types of events.This model, of parti- cular interest, describes among other things the input and departure processes of a finite-capacity Markovian node in a queueing network.
Even when P1 is Poisson, the analysis of run lengths is not entirely elementary.It requires careful attention to details, but once completed, it carries over easily to the MAP.In Section 2, we introduce the method for that case and obtain some explicit results.These are generalized to the MAP in Section 3 and the formal similarity of the derivations and results to the Poisson case is noted.Additional results for the MAP are derived in Section 4. Particular forms for the Poisson case are now simple special cases.Section 5 is devoted to a brief description of a computer-graphical exploration of the competition.Auxiliary results and the essential steps of lengthy derivations are found in the Appendix.

The Case of the Poisson Process
Let P1 be a Poisson process of rate .To visualize the competition between P1 and the sta- tionary grid, we superimpose the two point processes.We mark the endpoints of unit intervals during which r 1 or more Poisson events occur, and also the first Poisson events following strings of at least r 2 unit intervals devoid of Poisson events.These events mark the endpoints of 1-runs and 2-runs, respectively.
For a transparent analysis which immediately carries over to the MAP, we define an embedded Markov renewal process.That embedded process, which has an uncountably infinite state space, is conceptually easy.Its transition epochs are the marked endpoints.Those corres- ponding to the end of 1-runs (Poisson) are labeled with state 1.To those that are endpoints of 2- runs and occur at the time u during the corresponding unit interval, we associate the state (2, u).The bivariate sequence {(gn, rn)} of the states {Jn} at the successive marked epochs and of the sojourn times {rn} form a Markov renewal sequence on the state space E-{1} t.J {(2, u), with 0_<u<l}.
The case r 1 1 is somewhat special.we first assume that r I is at least two. of the Markov renewal process.
In order not to detract attention from the general case, Next, we construct the transition probability operator G The Transition Probability Operator While elementary, that construction involves mixed discrete-continuous probability mass functions and some absorption time distributions in an associated Markov chain.To prepare the way for the case of the MAP, we discuss the construction in detail.The operator G is specified by the conditional probabilities of transitions out of the state 1, and those of transitions out of a state of the type (2, u).Four cases are to be considered.
We use Gll(X) for the probability Gll(X) P{Jn + 1 1;rn + 1 -< X lJn 1}, for n _> 1 and x >_ 0. In a transition from state 1 to itself, that is, a l-run followed by another, the sojourn time is necessarily a positive integer ,.It equals if and only if the th interval of length one contains at least r I Poisson events, during all earlier unit intervals there are fewer than r 1 Poisson events, and there are no strings of r 2 1 consecutive intervals devoid of Poisson arrivals.The probability of that event is related to a simple (r 2 + 2)-state Markov chain with two absorbing states, whose transitions probability matrix P(rl, r2) is of the form here displayed for the representative value r 2 6.The non-zero elements are given by r 1 --1 a--e -X b(r1-1)-E e-)'-c(rl-1)-1-a-b(r1-1).
j=l At this time, we do not yet make the simplifications resulting from the special structure of the stochastic matrix P(rl,r2).It is preferable to emphasize its structural properties.We note that the transition probability matrix is of the form P(rl, r2) T(rl, r2) T o (r2) T,(rl) We see that, for x >_ 1, Gll(X) is the probability that, starting from the state 0, absorption into the state occurs before or at time Ix], where [. ]denotes the integer part.If we write /(r 2 -1), for the row vector of dimension r 2 -1 with components 1,0,..., 0, then [] Gll (x) E 7(r2-1)[T(rl'r2 )]j-1W'(rl )" 3--1 The sojourn times of transitions from state 1 to a state of the form (2, u') are necessarily of the form n + u', where n is an integer _> r 2 -1.To derive an explicit expression for the elementary conditional probability du,G12(u'; x) P{Jn + 1[(2, U'), (2, U' + du')]; 7 n + 1x Jn 1 }, we again consider the Markov chain with transition probability matrix P(rl,r2).At some time , that Markov chain must be absorbed in its state r 2. A run of r 2 events in P2 has now been completed.At some later time j, where j _ Ix-u'], the unit interval that contains at least one Poisson event must start and the first Poisson event must occur in the interval (u', u'+ du') count- ing from that epoch.In that case, there are j-unit intervals not containing Poisson events, in addition to the r 2 1 required for a 2-run.We see that [x .'1du'G12(u';x) E E 7(r2 1)[T(rl'r2)]-1w (r2)e-,k(j-)e-XU'du' 1 -e -' e .XU'Adu," Next, the transition probabilities starting from a state (2, u).First an observation that is obv- ious when P1 is Poisson, yet it is essential to the embedded process being a Markov renewal pro- cess.A corresponding property also holds when P1 is a MAP.It is then less obvious so we give an explicit argument in Lemma A-1 of the Appendix.Let the endpoint of the preceding 2-run fall at u in the unit interval during which it occurs.For the Poisson case, the number of additional events in P2 in that unit interval is then independent of the past and has Poisson distribution with parameter A(1-u).
If there are at least r 1 -1 Poisson events in that interval, the next event in the stationary unit grid is also the endpoint of a 1-run.In this case, we have a (2, u)---,l transition after a sojourn time of length 1 u.
If there are fewer than r I 1 Poisson events in that interval of length 1 u, then the sojourn times of the transitions, either to the state 1, or to some state (2, u'), are generated in exactly the same way as for the preceding two cases.
The expressions for the corresponding transition probabilities G21(u;x) and G22(u,u';x) are therefore as follows: since that mass-function has atoms only at points of the form 1-u + n, with n >_ 0. argument, 1 -e --X e u du'.
Laplace-Stieljes transforms of these mass-functions are evaluated in the Appendix.By setting s 0 in these transforms, we get the transition operator of the embedded discrete-time Markov process on the state space E. The Markov renewal process that concerns us here is a simple exam- ple of such processes on a general state space, see (inlar [1].For the present case, its analysis is entirely analogous to that of Markov renewal processes with finitely many states.The following statements are the salient results of that analysis.We note that (rl, r2) '(r 2 1)[I T(rl, r2) 1T'(rl), is the probability that the Markov chain with transition probability matrix P(rl, r2) is absorbed in the state . .It is also clear that 7(r 2 1)[I T(rl,r2) 1T (r 1 1 (rl, r2)   is the absorption probability into the state r 2.

tions and the normalization
The invariant probability measure on the state space E should satisfy the equa- 5) shows that 9()-Ke -u, for some constant K.
In the Appendix, it is shown that the fundamental mean of the Markov renewal process is given by This is an important quantity.Its inverse [E*(rl,r2)] -1 is the rate of transitions in the stationary version of the Markov renewal process.The ratio gl c(rl_ 1).I/1 (r 1, r2) E,(ri, r2) is the fraction of time spent in state 1 in the stationary version of the Markov renewal process.
Informally, it is the fraction of time in steady-state that the Poisson process is winning the competition to produce runs of its specified kind.Its inverse is the mean time between visits to State 1.
For completeness, we record the minor modifications needed when r 1 -1.There is no need to consider the absorption time in a Markov chain.The endpoints of unit intervals containing at least one Poisson event mark the ends of 1-runs, and the first Poisson event after a string of at least r 2-1 empty unit intervals marks the end of a 2-run.The Laplace-Stieltjes transforms of the transition operator for that case are:
Remark: The fact that, in all cases, l(rl, r2) equals c(r I -1) is to be expected.The visits to State 1 are the ends of 1-runs.These are precisely the endpoints of unit intervals during which at least r 1 Poisson events occur.The times between such points have a geometric distribution with mean [c(r 1 -1)]-1.
3. The Case of the Markovian Arrival Process We now generalize the preceding results to the case where P1 is a continuous-time Markovian arrival process with parameter matrices D O and D 1.We assume some familiarity with the basic properties of the MAP as discussed, e.g., in Lucantoni [2] and Neuts [3].The sum D of the para- meter matrices is an m m irreducible generator with stationary probability vector 0. The ma- trices P(n;t), n >_ 0, are the familiar matrices associated with the counting process.They satisfy the system of differential equations P'(0; t) P(0; t)D0, (9) P'(n; t) P(n; t)D o + P(n 1; t)D1, n >_ 1.
The sequence of matrices {P(n; t)} are a natural matrix generalization of the Poisson probabil- ity density.Methods for their numerical computation are discussed in Neuts and Li [5].To em- phasize the similarity to the Poisson probabilities and its role in the present analysis, we define the analogues of the quantities a, b(r 1 -1), and c(rl-1), but we use capital letters to emphasize that these objects are now m m matrices r 1 --1 A P(0; 1) exp(D0), B(r I 1) E P(j; 1), 3----1 C(r1-1)-E P(J;1) exp(D) A B(r1-1).

j=r I
The states of the continuous-parameter Markov chain with generator D are called the phases.
In modeling the competition for runs of lengths r I in P1 or r 2 in the reference process P2, the stationary unit grid, we construct a Markov renewal process analogous so that in Section 2. The only essential difference is that we must now also keep track of the phase at transition epochs.The state space is therefore E-{(1, j), 1 <_ j <_ m)t2 {(2, h, u), with 1 _< h _< m, 0 _< u < 1}.
For notational convenience, we denote the first set of states by 1, and a set of states with fixed u in the second by (2, u).
Lemma A-1 implies that the epochs labeled at the ends of runs, the future of the process is conditionally independent of its past, given the state in 1 or (2, u).The sequence of states and corresponding sojourn times therefore defines a Markov renewal process with the given states.We study the quantities associated with its stationary version as was done for the Poisson case.
The analogue of the stochastic matrix P(rl,r2) is now the matrix displayed here for the representative value r 2 6. consist of m states.The entries T'(rl) and T(r2) are therefore matrices of dimensions m(r2-1)m.By F(r2-1), we denote a matrix of dimensions mm(r2-1) of the form 1,0,...,0, where I is the identity matrix of order m.That matrix enables us to calculate simultaneously the absorption probabilities from each of the initial states (0, h), where h denotes the phase of the MAP.
The Transition Operator of the Markov Renewal Process For convenience, we give the expressions for the transition probabilities of the Markov renew- al process embedded at the ends of the runs immediately in transform forms.The derivation of the most complicated of these expressions, that in formula (12), is given in the Appendix; the others are similar but simpler.The analysis that follows applies when r 1 >_ 2 and r 2 _ 1.The simple special case where r 1 is treated at the end.
The matrix analogues of the transforms for the transition operator in the Poisson case are as follows: Gl(S e-SF(r2 1)[I-e-ST(rl,r2)]-1T'(rl).
Setting s-0 in the preceding formulas, we obtain the transition probability operator of the embedded Markov chain.We shall next derive its invariant probability measure explicitly.First, we introduce some matrices that are the analogues of the absorption probabilities in the finite Markov chain of Section 2.

K! Calculation of the Row Sum Means
As in the scalar case where P1 is the Poisson process, we need to calculate the row sum means of the semi-Markovian transition probability operator G.We proceed directly from the trans- forms in formulas (10-13).In each, we differentiate with respect to s, set s-0, and change sign.Next, we "sum" over the indices of the second state.For the transitions to the macro-state 1, that amounts to postmultiplication by a column vector e of ones.The results are column vectors mll and m21(u), 0 _< u < 1, all the dimension m.For transitions into a macro-state of the form (2, u'), we need to integrate between 0 and 1 with respect to u , followed by postmultiplication by e.These operations lead to column vectors m12 and m22(u), 0 _< u < 1, of dimension m.While the calculations are straightforward, they are substantially more complicated than in the Poisson case.They require thorough facility with the matrix formalism of the MAP.We state the final expressions in the form of a lemma.The salient calculational steps are given in the Appendix.

More Complex Descriptors
The descriptors of the MAP obtained so far, such as the mean sojourn times of the transitions between the states, are only moderately informative.Next, we study more complex descriptors based on a detailed analysis of the recurrence times of the states in the set 1. The time period between two successive visits to the set 1 will be called a 1-cycle.
The formal similarity of the general MAP to the special case of the Poisson process is now clear.Any special formulas for the Poisson case are obtained by using the scalar coefficient matrices D 1 -D O .Henceforth, the main discussion deals with the MAP only.
The expression for the fraction of time l(rl, r2) that the MAP wins the competition is en- tirely analogous to that of the Poisson case.It is given by gle Oc(r 1 1)e.
(37) l(rl,r2)--E*(rl,r2) This formula also has an intuitive explanation in terms of the mean of an arbitrary waiting time between unit intervals during which, in the MAP, at least r 1 events occur.

1-Cycles
The bivariate sequence of the phases at the starts of these recurrence times and their dura- tions is a Markov renewal sequence, which we shall call the process of 1-cycles.A typical 1-cycle consists either of a single transition from 1 to 1, or may be made up of a string of intervals that are the successive transitions from 1 to some state in (2, Ul) thence to (2, u2) and so on to a state in (2, un).From that state, there is a final interval in which the MAP wins and the end of the cycle is marked by a return to the set 1.
We shall study such descriptors as the number of intervals terminating in a state in (2, u), the fraction of time during the 1-cycle that the process P2 wins, and others.To that end, we define a number of matrices to enable us to write the transform formulas in concise forms that clearly bring out their structure.
For the transform of transition probability matrix of the process of 1-cycles, we use the notation E(,s).The element 2..,(,s) is the joint Laplace-Stieltjes transform of the conditional probability that the 1-cycle encs with the MAP in the phase j', that during the 1-cycle, the process P2 wins during a time at most x, the MAP wins during a time at most y, given that the 1-cycle starts with the MAP in the phase j.
To carry out the induction, it suffices to calculate the integral 1 J dvL*(n 1; , v)duG2(v u; ), 0 which, apart from matrices that can be factored out, reduces to the integral we have defined as the matrix B*(r I 1; , ).
The probability density (n;rl,r2) of the number of times the process P2 wins during an arbitrary 1-cycle, given by OC(r 1 1)n e (n;rl'r2) C(r1-1)e ( 50) for n >_ 0, is a useful, tractable descriptor.The quantity (n; rl, r2) is the probability that there are n gaps of length at least r 2 1 in a typical 1-cycle in the MAP.

Moments and Dependence
Clearly, the index n in En(,s is just a way of keeping track of another useful random variable, the number of 2-runs during the 1-cycle.The transform z) n--0 is immediately obtained in a closed form by noting that The calculation of various lower order moments from the transform is a matter of a routine (if involved) computation.It is possible, in principle, to express these quantities in terms of the coefficient matrices of the MAP.However, in doing so there are few analytic simplifications so that we prefer not to record the resulting formulas here.Without numerical algorithms, little insight can be gained from them, but the standard procedures for calculating these can be readily translated into algorithms for numerical computation.To compute the moment matrices, we differentiate twice in the matrix expression for (,s,z) and set s 0, and z 1 to obtain systems of linear equations.Their coefficients involve only vectors and matrices encountered earlier and are in a sufficiently modular form for numerical work.
In the cases where the 1-cycles can be expected to be long, say, because r is chosen fairly large, we can also expect that the behavior of the MAP during its successive 1-cycles will be only weakly dependent.The most readily available analytic tool to examine that dependence is the serial correlation of sequences such as the durations of successive 1-cycles and of the numbers of 2- runs during these.
Computationally useful expressions for these correlations are again routinely derived from the joint transform [OC(r 1 1)el-loC(rl-1)..(l,Sl,Zl)...(k, Sk, Zk)e of k successive 1-cycles.W(ll-known results for Markov renewal processes show that a readily available indication of weak correlation is how rapidly the successive powers of the stochastic mat- rix .(0, 0,1)-S(0, 0) converge to the limit matrix [OC(r 1)e]-leOC(rl-1).As that matrix is easily computed, that observation is useful in an initial numerical exploration of this issue.

Motivation and Applications
Before we undertook the present analytic study of the competition for runs we examined for some familiar MAPs by computer graphical displays.These suggested that this approach offers a conceptually simple way of bringing out important differences between MAPs.The classical des- criptors, based on second order properties, of these did not differ significantly.
For easily simulated stationary point processes, such as the MAP, the data needed to high- light the competition are easily collected.It is sufficient to attach the value of a counter to each event in the given process and in the reference process.That counter, which can conveniently take positive integer values for the first and negative ones for the second, keeps track of how many uninterrupted events of each type have currently occurred.For these counts, the 1-runs and 2-runs corresponding to various values of r 1 and r 2 are easily identified.The effect of these parameters on a long simulation record can be examined without having to repeat the simulation and without making several passes through a long record of data.
To simplify interpretation of the results, we generated streams of several thousand events for various stationary MAPs of unit rate.The streams of interest and the reference process (the grid, in the present case) all had a common rate.We fixed r I at some integer value of modest size, such as eight or ten.For each event stream, we graphically recorded the ends of the successive 1- runs and 2-runs by printing a symbol 1 or 2. That was done for several increasing values of r 2.
The resulting graphs and elementary summary statistics provided clear pictures that were remarkably consistent for different simulations with the same parameter matrices.In contrast, differences between moderately distinct MAPs stood out vividly.The successive 1-cycles typically involved several hundreds of events and took on a small number of typical visual forms.The prevalence of gaps (as measured by the occurrence of 2-runs during 1-cycles) stood out clearly and, as is to be expected, decreased with increasing values of r 2. These, and the occurrence of bursts in affecting the frequency of 1-runs, clearly distinguished between the MAPs.
We concluded that the competition for runs brings out salient local properties of the realizations of these point processes, features that are missed by the averaging inherent in descriptors based on second order properties or on spectral methods.While that conclusion is based on computer experiments, we are confident that numerical computations of the descriptors in Section 4 will substantiate these findings.The present results, and others to be included in sub- sequent articles, should serve in a formal quantification of these experimental perceptions.
The results in Section 4 show that analytic use of the competition for runs as a tool in the in- vestigation of MAPs requires considerable numerical work.The transform formulas obtained there clearly indicate which embedded Markov renewal processes need to be studied in detail.
However, their complex forms show that this is best done by modular numerical algorithms than by explicit formulas.Computer implementation remains to be done, but the resulting, well- computable descriptors should prove very useful.
The input and departure processes of a finite-capacity Markovian node in a queueing network, to which we alluded in Section 1, offer examples of two dependent processes for which the competition for runs is worth examining.The features of the competition bring out how the number of jobs within the node "pulsates" over time.Extensive numerical computations for small nodes, based on the MAP with two event types, should serve to quantify these features for benchmark examples.For less structured stationary point processes, it appears that only detailed simulation studies can indicate whether the qualitative findings for MAPs have wider applicability.
In the postmultiplication by e, we note that, since D-D 0+D1 is a generator, (-D0)-ID1 e-e.That results in the cancellation of several terms and formula (31) follows after routine simplifications.The verification of (33) proceeds along the same lines.