**Abstract** - We start with the premise, and provide evidence that it is valid, that a Markov-modulated Poisson process (MMPP) is a
good model for Internet traffic at the packet/byte level. We present an algorithm to estimate the parameters and size of a discrete MMPP (D-MMPP) from a data trace. This algorithm requires only two passes through the data. In tandem-network queueing models, the input to a downstream queue is the output from an upstream queue, so the arrival rate is limited by the rate of the upstream queue. We show how to modify the MMPP describing the arrivals to the upstream queue to approximate this effect. To extend this idea to networks that are not tandem, we show how to approximate the superposition of MMPPs without encountering the state-space explosion that occurs in exact computations. Numerical examples that demonstrate the accuracy of these methods are given. We also present a method to convert our estimated D-MMPP to a continuous-time MMPP, which is used as the arrival process in a matrix-analytic queueing model.

**Index Terms** - Hidden Markov model, Markov-modulated Poisson process (MMPP), matrix-analytic queueing model, superposition, tandem queues.

### I. INTRODUCTION

Two of the many many applications of queueing models in communications networks are sizing links in transport networks and buffers in routers. A fundamental part of a queueing model is the arrival process. A Markov-modulated Poisson process (MMPP) is an attractive model for describing backbone packet traffic. This paper describes a computationally simple method of fitting an MMPP to trace data. Suppose several traffic streams, each described by an MMPP, are going to be combined. This occurs when access links are combined at a provider edge router and when several links are replaced by a
single higher-speed link. This paper provides an algorithm for combining the MMPPs for the original streams into a tractable MMPP that describes the superposition of these streams.
Once it was recognized that (at the packet or byte level) IP traffic is not described well by a Poisson process, it was natural to try and describe it by an MMPP (see, for example, [1]-[3]). The MMPP has enough flexibility to describe a wide variety of data, and its physical interpretation seems to describe rate fluctuations in many situations. For example, there are several methods of fitting an MMPP to capture the behavior over several time scales (see, e.g., [4] and [5]). Therefore, the self-similar properties described in [6] can be matched over the time scales of interest. Moreover, the time-varying Poisson properties of the MMPP is consistent with the statistical findings of [7] for data collected on highly utilized links. In addition, there is a wide body of knowledge and software available to make computable performance models [8], [9]. A challenge in describing a time
series by an MMPP is to estimate the parameters of the model accurately while keeping the number of states small enough to make the performance models tractable. This paper addresses that challenge by proposing an algorithm to estimate MMPP parameters from count data.

An additional challenge is to limit state-space explosion in the description of a superposition of MMPPs. The superposition of two MMPPs is an MMPP whose transition matrix is the Kronecker sum of the transition matrices of the constituent MMPPs. The dimension of the superposition of *κ* MMPPs with dimensions *η*_{1}, η_{2}...,η_{k} is the product of these dimensions, so performance models where the arrival process is the superposition of MMPPs become intractable *κ* for as small as two when the original MMPPs have about 20 states. We use the idea behind the MMPP estimation algorithm to approximate the superposition MMPP.We find that the number of states grows slowly with the number of MMPPs superposed.

A third contribution of this paper is motivated by the following scenario. Users generate packet arrivals in the form of anMMPP. The arrivals pass through a constant-rate device, such as an access line (so many bits per second) or a router (so many packets per second). When the user rate exceeds the rate of the device, assume that the excess is either buffered or retransmitted after a short delay. Then the output from the device is a modification of the original MMPP with a peak rate equal to the rate of the device. We present a method to capture this effect and develop
equations to describe this rate-limited process.

With the discovery that much Internet traffic exhibits longrange dependence (LRD), one might wonder why we propose to use an MMPP that, by definition, is short-range dependent (SRD), as a model of IP traffic. As discussed in [10], although SRD traffic has an exponential tail, the exponential part of the tail might not become dominant until far beyond the region that is of interest for performance modeling. For example, in [11],
realistic examples are given where an MMPP has a "heavy tail" until beyond the point where the probability of blocking is less than 10^{-40}. Beyond that, the exponential component of the tail is dominant but since the model is applied to an ATM traffic example, we are only interested in regions where blocking is no less than 10^{-9}. There has also been some recent work [4], [5] on fitting an MMPP to IP traffic such that the LRD behavior is matched over several orders of magnitude.

Also, the state-of-the-art in computing performance measures with MMPP models is far beyond that of current LRD modeling. For example, a recent book on the state-of-the-art of LRD modeling was edited by K. Park and W.Willinger. In their overview paper, [6], they summarize the state-of-the-art as illustrated with the following three quotes: "A major weakness of many of the [LRD] queueing based results is that they are asymptotic, in one
form or another," "A further drawback of current performance results is that they concentrate on first-order performance measures that relate to (long-term) packet loss rate but less so on second order measures - for example variance of packet loss or delay, generically referred to as jitter - which are of importance in multimedia communication" and finally, "Even less is known about transient performance measures, which are more relevant
in practice when convergence to long-term steady-state behavior is too slow to be of much value for engineering purposes." We note that results for the MMPP/G/1 queue are exact, that results exist for moments as well as distributions and that exact results for the transient distributions are also available [8], [12].

This paper is organized as follows. In Section II, we present a method to fit a discrete MMPP (D-MMPP) to an IP data trace. In Section III, we present a method to approximate the effect of limited access line speeds on the traffic. Building upon both of these ideas, we present in Section IV a method for approximating the superposition of D-MMPPs with a D-MMPP having a much smaller state space. This will allow an analytic treatment of the multiplexing of several sources (some of which might represent outputs from previous nodes). Section V contains several numerical examples demonstrating the accuracy of these methods both in terms of approximating the arrival process as well as predicting the performance seen by feeding the process into a queue. Finally, we conclude the paper in Section VI.

### II. FITTING AN IP TRACE

Heffes [13] made the first documented proposal to fit a traffic stream by an MMPP. He used two phases for the Markov chain, gave an estimation procedure for the parameters of the MMPP from count data, and solved several queueing models with exponential service times. Meier-Hellstern [14] provided the first estimation algorithm based on interarrival-time measurements. It has proven difficult to obtain estimators for MMPPs with more
than a few states because of the computational burden. Moreover, traffic data frequently consists of counts of events during fixed length intervals, such as packets per second or bytes per 100 ms, so model fitting has to be done with this type of data.
Hidden Markov models take counts as the fundamental random variables, so they are alternatives to MMPP models. A *hidden Markov model* (HMM) is a discrete-time process in which parameters that control the distribution of the number of counts depend on the state of a (hidden) Markov chain [15].
In particular, when the number of counts has a Poisson distribution, we obtain a *Poisson-hidden Markov model* (PHMM). The connection between the MMPP and the PHMM is not straightforward. An MMPP is a point process and the fundamental random variables are the arrival epochs. The phases of
an MMPP are governed by a *continuous-time* Markov chain. Given the current phase, the future of the MMPP is conditionally independent of the past (i.e., Markov); this is what makes the MMPP so useful in continuous-time queueing models.

Suppose we define a point process where the phases are governed by a *discrete-time* Markov chain and the rates depend on the phases exactly as in an MMPP. The phase changes can occur at integer multiples of some fundamental time unit; call these times *phase-transition* epochs. This process is not Markovian because the times between phase changes do not have the memoryless property; it is regenerative (at any phase-transition epoch). The essential physics are the same as the MMPP. We propose calling this process a *discrete* MMPP (D-MMPP). A PHMM is a counting process that counts the number of events between adjacent phase-transition epochs of a D-MMPP. It is a special case of a batch Markovian arrival process (BMAP); Lucantoni [8] describes the BMAP and some related processes.

With the counts as data *and the number of states in the Markov chain given*, the log-likelihood function of a PHMM can be written in a computationally useful form. Moreover, it is claimed that extant optimization algorithms can be used to obtain maximum likelihood estimates [15]. The sticking point in this result is that the number of Markov chain states has to be specified in advance, and that may not be easy to do. The claim
that maximizing the likelihood function can be done needs to be checked for time series with hundreds of thousands of observations and up to forty states in the Markov chain, which may be required for IP traffic counts. Doing so is beyond the scope of this paper.

Previous approaches to estimating MMPP parameters when there are more than two states use the maximum likelihood estimator (MLE) method [16]. The MLE is also used to estimate parameters of the PHMM [15]. MLEs are computationally expensive because a nonlinear function is minimized; each function evaluation has computational effort that is linear in the length of the trace that is being fitted, and quadratic in the number of states in the Markov chain underlying the MMPP. This may not be tractable for 10 000 observations and 20 states, which is what some of our traffic data yields. (One count every second for three hours produces 10 800 counts). Moreover, the MLE method requires that the number of states of the Markov chain is picked before the calculations are done. Algorithm LAMBDA in Section II-B1 can be used to do this.

The interrupted Poisson process (IPP) is the simplest nondegenerate special case of an MMPP. There are two states, and the arrival rate is zero in one of them. Andersson and Ryden [16] developed a maximum likelihood estimation algorithm for MMPPs based on the superposition of IPPs. Their algorithm requires that the arrival times are given in addition to the counts. We want to avoid needing that data. Lee *et al.* [17] describe the properties of a generalized IPP in which the times between phase transitions need not be exponential. They focus on the autocorrelations in the times between arrival epochs. Since only one positive Poisson rate is used, we do not think this model will fit the data we use (see Fig. 1.)

**A. D-MMPP Model**

The D-MMPP is defined on a discrete Markov chain with. . . **(continue to PDF version)**