3.2.4.1  Long-range dependence

The process X is long-range dependent if, for every fixed T_i, the autocorrelation function is not summable:

3.27 $$\sum _{k=0}^{\infty }\text{Cov}(k)=\infty$$

A better definition is that a process with LRD has a covariance function with hyperbolic decay:

3.28 $$\text{Cov}(X_n,X_{n+k})\sim |k{|}^{-\beta }$$

when |k| → ∞ and 0 < β < 1 (the definition of parameter β is given by Equation 3.23.

On the contrary, for the SRD processes, the covariance function has an exponential decay:

3.29 $$\text{Cov}(X_n,X_{n+k})\sim a^{\left|k\right|}$$

when |k| → ∞ and 0 < a < 1. Therefore, a SRD process has a summable autocorrelation function

3.30 $$\sum _{k=0}^{\infty }\text{Cov}(k)=\text{finite}<\infty$$

We can conclude that Cov(k) becomes small for big values of k, while their cumulative effect is significantly more important in the case of the LRD processes as the one shown in Figure 3.1.

3.2.4.2  Slowly decaying dispersion

As was already mentioned, the dispersion of the samples average decays slower in the case of self-similar processes for m big enough:

Figure 3.1   Autocorrelation structures for SRD processes (a) and for LRD processes (b).

3.31 $$\text{Var}[X^{(m)}]\sim m^{-\beta }$$

On the other hand, for SRD processes β= 1 and

3.32 $$\text{Var}[X^{(m)}]\sim m^{-1}$$

The property of slowly decaying dispersion can be easy shown by the graphic plot of Var[X^(m)] as function of m on a logarithmic axis (dispersion–time diagram). A straight line with negative slope <1 (on a large range of values of m) points to a slowly decaying dispersion. The same property can be expressed by the dispersion index of counting F(T) (given in Equation 3.7).

The equation

3.33 $$\text{Var}[X^{(m)}]=\frac{E[N(m{T}_i)]F[m{T}_i]}{m^2}$$

shows the equivalence between a dispersion–time diagram, which indicates a slowly decaying dispersion with parameter β and a plot of the dispersion index with slope 1 −β.

3.2.4.3  Heavy-tailed distributions

If the abovementioned properties (covariance and autocorrelation) refer the scaling properties of time-dependent statistics, the attributes of the heavy-tailed distributions are focussed on the marginal amplitudes distribution of X_n, as defined in Equation 3.1. Packets arrivals or ON/OFF intervals are examples of such processes.

A random variable X is called heavy-tailed if the complementary cumulative distribution function (CCDF) decays as a power law:

3.34 $$\mathrm{Pr}P\{X\ge x\}\sim x^{-\alpha }$$

when x → ∞ for α > 0. For 0 < α ≤ 1 all the moments of X are infinite and in a wide sense, the nth moment is infinite for n ≥ α. The most known type of heavy-tailed distribution is the Pareto distribution. Let us also mention that although the property of heavy-tailed distribution is not a necessary condition for self-similarity, the self-similar nature of many traffic data results from heavy-tailed distributions (the so called Joseph effect, which illustrates the idea that movements in a time series tend to be part of larger trends and cycles more often than they are completely random; the Joseph effect is quantified by the Hurst component, where movements fall between a Hurst range of 0 and 1).

3.2.4.4  1/f noise

This property is an aspect of the LRD in the frequency domain, more precisely the evidence of the fact that the density of the power spectrum (PSD—defined in Equation 3.8) of the LRD processes diverges at small frequencies, as a contrary in comparison with the SRD processes, which have a smooth PSD at small frequencies, that is,

3.35 $$S(\omega )\sim \frac{1}{|\omega {|}^{\gamma }}$$

when |ω| →0 and 0 <γ< 1. The relation between the parameters γ and β is

3.36 $$\gamma =1-\beta =2H-1$$

For SRD processes, which have a finite PSD when |ω| →0, (i.e., when γ= 0), the values of the autocorrelation function R(k) decays fast for big values of k, reaching a finite sum.

3.2.4.5  Fractal dimension

The fractal dimension d of an object is defined by

3.37 $$d=\frac{\ln N}{\ln (1/\eta )}$$

where N is the number of self-similar pieces that cover an d-dimensional object, and η is the dimension of the scale unit. For SRD processes, d = 1, while for LRD processes, d has fractioned dimension of form d = xx.yy.….

3.2.4.6  SRD versus LRD: Conclusion

For LRD processes, the following characteristics are representative:

• $E[X_n^m]$
  is close to the second-order moment when m → ∞.
• Var[X^m] approaches asymptotically Var[X]/m when m → ∞.
• ${\sum }_{k=0}^{\infty }\text{Cov}(X_n,X_{n+k})$
  is convergent.
• The spectrum S(ω) is finite for ω= 0.

For LRD processes, the following characteristics are representative:

• $E[X_n^m]$
  differs the second-order moment when m → ∞.
• Var[X^m] approaches asymptotically m^−β when m → ∞.
• ${\sum }_{k=0}^{\infty }\text{Cov}(X_n,X_{n+k})$
  is divergent.
• The spectrum S(ω) is singular for ω= 0.

3.3.1.1  Pareto distribution

Pareto distribution is a power-law function parameterized by two parameters: the shape parameter α (named also Pareto index) and the scale parameter β (named also cutoff parameter because it represents the minimum value of the random variable X). The cumulative distribution function (CDF) Pareto is given by

3.38 $$F(x)=1-{\left(\frac{\beta }{x}\right)}^{\alpha }$$

and the probability density function for x > β and α > 0

3.39 $$f(x)=\frac{\alpha }{\beta }{\left(\frac{\beta }{x}\right)}^{\alpha +1}$$

for x ≤β

3.40 $$f(x)=F(x)=0$$

The mean (the expected value) is

3.41 $$E[X]=\frac{\alpha \beta }{\alpha -1}\text{with}\alpha >0$$

The dispersion is

3.42 $${\sigma }^2=\frac{\alpha }{{(\alpha -1)}^2(\alpha -2)}$$

with α> 2 and β= 1. The parameter α determines the mean and the dispersion as follows:

• For α≤ 1, the distribution has infinite mean.
• For 1 ≤α≤ 2, the distribution has infinite mean and infinite dispersion.
• For α≤ 2, the distribution has infinite dispersion.

The relation between parameter α and Hurst parameter is

3.43 $$H=\frac{3-\alpha }{2}$$

3.3.1.2  Log-normal distribution

A log-normal distribution is a probability distribution of a random variable whose logarithm is normally distributed. If the stochastic variable Y = logX is normally distributed, then the variable X is log-normally distributed. The distribution has two parameters, the mean μ of ln(x), called also location parameter (μ> 0) and the standard deviation σ of ln(x), called also the scaling parameter (σ > 0), for 0 ≤ x ≤ ∞. The cumulative distribution function and the probability density function are

3.44 $$F(x)=\frac{1}{2}\left[1+\text{erf}\left(\frac{\ln x-\mu }{\alpha \sqrt{2}}\right)\right]$$

3.45 $$f(x)=\frac{1}{\alpha x\sqrt{2\pi }}{\text{e}}^{-{(\ln x-\mu )}^{2}2{\sigma }^2}$$

The error function erf(y) is defined as

3.46 $$\text{erf}(y)=\frac{2}{\sqrt{\pi }}\underset{0}{\overset{y}{\int }}{\text{e}}^{-t^2}\text{d}t$$

The mean and the dispersion of the log-normal distribution are

3.47 $$E[X]={\text{e}}^{\mu +({\sigma }^2/2)}$$

3.48 $$\text{Var}[X]={\text{e}}^{2\mu +{\sigma }^2}({\text{e}}^{{\sigma }^2}-1)$$

The log-normal distributions are useful in many applications such as regressive modeling and the analysis of experimental models, especially because it can model errors produced by many factors.

3.3.2.1  ON–OFF models

This approach was first suggested by Mandelbrot, for economical maters, based on renewal processes and later was applied by W. Willinger, M. Taqqu, R. Sherman, and D. Wilson for network communication traffic (Willinger et al. 1997).

The aggregate traffic is generated by multiplexing several independent ON–OFF sources, with each source changing these states in a given periodicity. The generated traffic is (asymptotically) self-similar due to the superposition of several renewal processes, where the renewal values are only 0 and 1, and the interval between renewal is heavy-tailed (e.g., a Pareto distribution with 1 <α< 2). This superposition presents the syndrome of the infinite dispersion (the so-called Noah effect) and leads to closing of the fractional Brownian motion (fBm) in self-similar aggregate traffic when the number of sources is very high (the so-called Joseph effect), with the Hurst parameter being

3.49 $$H=\frac{3-\alpha }{2}$$

The ON–OFF model is simple and therefore very attractive for the simulation of Ethernet traffic, considering each station connected in an Ethernet network as a renewal process, staying in state ON when communicating and in the state OFF when not communicating. So Ethernet traffic is the result of aggregation of several renewal processes.

3.3.2.2  Fractional Brownian motion

The concept of fractional Brownian motion (fBm) was introduced by Mandelbrot and Van Ness in 1968 like a generalization of the Brownian motion (Bm). The Bm(t) process has the property that any random vector at different times has a joint distribution identical to that of a rescaled and normalized version of the random vector, which is a unique continuous Gaussian process with this property and with stationary increments. The (standard) fBm family of processes B_H(t) is indexed by a single parameter H ∈ [0, 1], which controls all of its properties. Larger values of H correspond to stronger large-scale correlations, which make the sample paths appear more and more smooth. In contrast, B_1/2 is the Brownian motion without memory, which in discrete time is just a random walk. The increments of fBm with H∈ [0, 1] define the well-known fractional Gaussian noise (fGn) family of stationary processes. For Brownian motion, these increments are i.i.d., which is clearly SRD. fGn processes with H∈ [0, 1/2) are also SRD, though qualitatively different from white noise, whereas fGn processes with H∈ (1/2, 1] are LRD. Of course, this appearance can be modified by allowing for two additional parameters, namely mean and variance, which allows capturing a wide range of “smoothness” scenarios, from “highly bursty” to “very smooth.” For analytical studies, the canonical process for modeling LRD traffic rate processes is fGn with H∈ (1/2, 1]. When describing the observed self-similar behavior of measured traffic rate, we can be more flexible. The standard model is a self-similar process in the asymptotic sense or, equivalently, an LRD process.

The following definitions are specific for a fBm process:

The expected value of the process is zero

3.50 $$E[B(t-\tau )-B(t)=0]$$

The dispersion is proportional to the time difference

3.51 $$\text{Var}[B(t+\tau )-B(t)]={\sigma }^2\left|\tau \right|$$

The fGn process B_H(t) is defined as the moving average of dB(t), where the past increments of B(t) are mediated by the kernel (t−τ)^H − (1/2) and can be expressed by the integral fractional Weyl integral of B(t):

3.52 $$B_H(t)=\frac{1}{\Gamma (H+(1/2))}\underset{-\infty }{\overset{0}{\int }}{(t-\tau )}^{H-(1/2)}-{(-\tau )}^{H-(1/2)}\text{d}B(\tau )+\underset{0}{\overset{t}{\int }}{(t-\tau )}^{H-(1/2)}\text{d}B(\tau )$$

where Γ(x) is the Gamma function. The value of a fBmprocess at moment t is

3.53 $$B_H(t)=X{t}^H$$

where Xis a random variable normally distributed with zero mean and dispersion equal to 1, where due to the self-similarity

3.54 $$B_H(at)=a^H{B}_H(t)$$

For H = 0.5, one obtains the ordinary Brownian motion (Bm).

To conclude, let consider the main statistical properties of a fBm process

• Mean: E[B_H(t)] = 0
• Dispersion: Var[B_H(t)] = Var[Xt^H] = t^2H
• Autocorrelation: R_BH(t, τ) = E[B_H(t)B_H(τ)] = 1/2(t^2H +τ^2H − |t −τ|^2H )
• Stationary increments: Var[B_H(t) − B_H(τ)] = |t −τ|²

3.3.2.3  Fractional Gaussian noise

The stationary process of increment X_H of fBm type is known as fractional Gaussian noise (fGn):

3.55 $$X_H=(X_H(t)=B_H(t+1)-B_H(t):t\ge 0)$$

It is a process with Gaussian distribution, with mean μ and dispersion σ². The dispersion of the mean of the square of increments is proportional to the time difference, and the autocorrelation function is

3.56 $$R(t)=\frac{1}{2}(|t+1{|}^{2H}-|t{|}^{2H}+|t+1{|}^{2H})$$

In the limit case when t → ∞

3.57 $$R(t)\sim H(2H-1)|t{|}^{2h-2}$$

A fGn process is exactly self-similar when 0.5 < H < 1. When μ_X = 0, fGn becomes Gaussian time-continuous process, $B_H={\{B_H(t)\}}_{t=0}^{\infty }$

3.3.2.4  Fractional ARIMA processes (FARIMA)

An integrated autoregressive moving-average (p, d, q) series, denoted ARIMA (p, d, q), is closely related to ARMA (p, q). ARIMA models are more general as ARMA models and include some nonstationary series (Box and Jenkins 1976).

FARIMA processes are the natural generalizations of standard ARIMA (p, d, q) processes when the degree of differencing d is allowed to take nonintegral values. The mathematical aspects were discussed in detail in Section 2.5.5.

Let remember that for 0 < d < 0.5, FARIMA (p, d, q) process is asymptotically self-similar with H = d + 0.5. An important advantage of this model is the possibility to capture both short- and long-range-dependent correlation structures in time series modeling. The main drawback is the increased complexity, but it can be reduced using the simplified FARIMA model (0, d, 0).

Consequently, compared with ARMA and fGn processes, the FARIMA (p, d, q) processes are flexible and parsimonious with regard to the simultaneous modeling of the long- and short-range-dependent behavior of a time series.

3.3.2.5  Chaotic deterministic maps

A chaotic map f:V →V can be defined by the following three criteria (Erramilli et al. 2002):

1. It shows sensitive dependence on initial conditions (SIC), which means that typical trajectories starting from arbitrarily close initial values nevertheless diverge at an exponential rate. SIC is the basis of the “deterministic” property of chaotic systems, because errors in the estimates of initial states become amplified and prevent accurate prediction of the trajectory’s evolution.
2. It is topologically transitive or strongly mixing, which implies that it is irreducible and that it cannot be decomposed into subsets that remain disjoint under repeated action of the map.
3. Periodic points are dense in V.

From a modeling perspective, even low-order chaotic maps can capture complex dynamic behavior. It is possible to model traffic sources that generate fluctuations over a wide range of timescales using chaotic maps that require a small number of parameters.

If one considers a chaotic map defined as X_{n + 1} = f(X) and two trends (evolutions) with almost identical conditions X₀ and X₀ + ε, then SIC can be expressed by

3.58 $$|f^N(X_0+ϵ)-f^N(x_0)|=ϵ{\text{e}}^{N\lambda (x_0)}$$

where f^N(X₀) refers to a N time-iterated map, and the parameter λ(x₀) describes the exponential divergence (the so-called Liapunov exponent). For a chaotic map, λ(x₀) should be positive for “almost all” X₀. In other words, the limits in the accuracy of the specification of the initial conditions lead to an exponential growth of the incertitude, making the traffic behavior on long-term unpredictable. Another important characteristic is the convergence of the trajectories (in the phase space) to a certain object named strange attractor, which is typical for a fractal structure.

In order to demonstrate how we can associate packet source activity with the evolution of a chaotic map, let consider the one-dimensional map in Figure 3.2 where the state variable x evolves according to f₁(x) in the interval [0, d) and according to f₂(x) in the interval [d, 1)].

Figure 3.2   Evolutions in a chaotic map.

One can model an on–off traffic source by stipulating that the source is active and generates one or more packets whenever x is in the interval [d, 1)], and is idle otherwise. By careful selection of the functions f₁(x) and f₂(x), one can generate a range of on–off source behavior. The traffic generated in this manner can be related to fundamental properties of the map. As the map evolves over time, the probability that is in a given neighborhood is given by its invariant density ρ(x), which is the solution of the Frobenius–Perron equation:

3.59 $$\rho (x)=\underset{0}{\overset{1}{\int }}\delta (x-f(z))\rho (z)\text{d}z$$

If the source is assumed to generate a single packet with every iteration in the on-state (d ≤ x ≤ 1), the average traffic rate is simply given by

3.60 $$\lambda =\underset{0}{\overset{1}{\int }}\rho (x)\text{d}x$$

Another key source characteristic is the “sojourn” or run times spent in the on- and off-states. Given the deterministic nature of the mapping, the sojourn time in any given state is solely determined by the reinjection or initial point x₀ at which the map reenters that state. For example, let the map reenter the on-state at x₀ ≤ d. The sojourn time there is then the number of iterations it takes for the map to leave the interval [d, 1]:

3.61 $$\begin{array}{cc}\hfill f_2^{k-1}(x_0)\le d,\hfill & \hfill f_2^k(x_0)<d\hfill \end{array}$$

One can then derive the distributions of the on- and off-periods from the reinjection density ψ(x) (where the probability that the map reenters a state in the neighborhood [x₀, x₀ + dx] is given by ψ[x₀]dx) and Equation 3.61. The reinjection density can be derived from the invariant density

3.62 $$\psi (x)=\rho (f^{-1}(x))\frac{\text{d}{f}^{-1}(x)}{\text{d}x}$$

One can use these relations to demonstrate the following.

1. The sojourn times in the two states are geometrically distributed if f₁(x) and f₂(x) are linear:
   3.63 $$\begin{array}{cccc}\hfill f_1(x)=\frac{x}{d},\hfill & \hfill 0\le x<d;\hfill & \hfill f_2(x)=\frac{1-x}{1-d},\hfill & \hfill d\le x\le 1\hfill \end{array}$$

   
   This follows from the fact that the invariant density in this case is uniform, and substituting for the linear mappings in Equation 3.60.
2. In order to match the heavy-tailed sojourn time distribution behavior observed in measurement studies, one can choose f₁(x) such that as x → 0, f₁(x)~ x + cx, where 1.5 < m< 2. Note that this function evolves very slowly in the neighborhood of 0 (a fixed point). This results in sojourn times in the off-state that are characterized by distributions with infinite variance.

One can similarly generate extended on-times by choosing f₂(x) so that it behaves in a similar manner in the neighborhood of 1 (another fixed point). In this way, one can match the heavy-tailed sojourn time distribution observed in measurement studies.

Thus, chaotic maps can capture many of the features observed with actual packet traffic data, with a relatively small number of parameters. While the qualitative and quantitative modeling of actual traffic using chaotic maps is itself the subject of ongoing research, ultimately the usefulness of this approach will be determined by the feasibility of solving for the 1-D and 2-D invariant densities. We have demonstrated the potential of chaotic maps as models of packet traffic. Given that packet traffic is very irregular and bursty, it is nevertheless possible to construct simple, low-order nonlinear models that capture much of the complexity. While chaotic maps allow concise descriptions of complex packet traffic phenomena, there are considerable analytical difficulties in their application.

As we indicated earlier, a key aspect of this modeling approach is to devise chaotic mappings f(x) such that the generated arrival process matches measured actual traffic characteristics, while maintaining analytical tractability. Other source interpretations relating the state variable to packet generation are possible and may be more suitable for some applications.

Performance analysis, that is, analysis of queuing systems in which the randomness in the arrival and/or service characteristics is modeled using chaotic maps, particularly systems without a conventional queuing analog requires additional work on the modeling and representation of systems of practical interest using chaotic maps, as well as investigation of numerical and analytical techniques to calculate performance measures of interest, including transient analysis.

3.3.3.1  ON–OFF method

This method (also called aggregation method) was described in Section 3.3.2. It consists in the aggregation and superposition of renewal rewards process (ON/OFF), in which activity (ON) and inactivity (OFF) periods follow a heavy-tailed PDF. We also know that there are some drawbacks associated with the deployment of this technique in the network simulation procedures, such as pitfalls in choosing timescales of interest and number of samples. Despite these issues, this approach could allow an immediate use of widespread network simulation tools, such as Network Simulator 2 – ns2 or the software family from OPNET since there is no need to extend their libraries to support such analytical models.

The main advantage is the possibility to implement the algorithm for parallel processing. The main drawback of the aggregation technique is its large computational requirements per sample produced (because of the large number of aggregated sources). This is the main reason why parallelism is introduced. However, even though the computational requirements are large, its actual computational complexity is O(n) for n produced samples. All the other self-similar traffic generation techniques exhibit poorer complexity than O(n) and they eventually result in a slowdown of the simulation as the simulated time interval increases. As an example, Hosking’s method requires O(n²) computations to generate n numbers (Jeong et al. 1998).

To illustrate the above statements, let consider the following simulation model. Cell arrivals will be represented by Run Length Encoded (RLE) tuples. An RLE tuple t_i includes two attributes, s(t_i), the state of the tuple, and d(t_i), the duration of the tuple. The two attributes represent the discrete time duration d(t_i) over which the state s(t_i) stays the same. The state is either an indication of whether the source is in the ON or OFF state (e.g., 0 for OFF and 1 for ON in a strictly alternating fashion), or the aggregate number of sources active (in the ON state) for the specified duration, that is, for N sources, s(t_i) ∈ {0, 1, …, N}. Thus, a sequence of t_i’s is sufficient for representing the arrival process from an ON/OFF source or from any arbitrary superposition of such sources. The benefit of such representation is that the activity of the source over several time slots can be encoded as a single RLE tuple. In summary, the algorithm proceeds by generating a large number, N, of individual source traces in RLE form. The utilization of each one of these sources is set to U/N, such that the aggregation of their N results in the desired link utilization to U. Each logical process (LP) of the simulation merges and generates the combined arrival trace for a separate nonoverlapping segment of time, which we call a slice. Thus, each LP is responsible for the generation of the self-similar traffic trace in the form of RLE tuples over a separate segment (slice) of time.

In logical terms, the concatenation of the slices produced by each LP in the proper time succession is the desired self-similar process. The LPs continue looping generating a different slice each time. The LP performs the generation of the self-similar traffic trace by going through the following three steps at each simulated time slice:

1. It generates the merge of the RLE tuple traces of the N individual sources.
2. It aggregates the merged traffic into a link speed equal to the desired access link speed.
3. It corrects the produced RLE trace by incorporating any residual cell counts. The generation of the individual RLE source traces is also performed in parallel. Each LP generates all the slices of a subset of sources that will be necessary to a number of P LPs during the generation of the current slices. That is, if P LPs are participating in the simulation, each LP generates the RLE tuples of P subsequent slices for the sources that it has been assigned to generate. It then sends their P −1 to the other LPs for each source it simulates. The individual ON/OFF sources are parameterized accordingly to fit the desired self-similar traffic. Namely:
   • The shape value, α, of the Pareto distribution used for the ON period is set according to H = (3 −α)/2, where H is the desired Hurst value.
   • Since N ON/OFF sources are aggregated, the per-source utilization of U/N is determined by the ratio of the ON and OFF periods of the individual processes. That is, the average OFF period E[OFF] is set to E[OFF] = E[ON].
   • The average ON period E[ON] is set to E[ON] = B, the average burst length, which can be derived from traffic measurements. E[ON] does not have any impact on the self-similarity, and it can be considered a free variable.

3.3.3.2  Norros method

An exactly self-similar can be generated using a fBm Z(t) process

3.64 $$\begin{array}{cc}\hfill A(t)=mt+\sqrt{am}Z(t)\hfill & \hfill \text{for}-\infty <t<\infty \hfill \end{array}$$

A(t) represents the cumulative arrivals process, that is, the total traffic that arrives until the time t. Z(t) is a normalized fBm process with 0.5 < H < 1.0, where m is the average rate, and a is the shape parameter (the ratio between dispersion and mean in the time unit) (Norros 1995). We call this an (m, a, H) fBm. When this arrival process flows into a buffer with a deterministic service, the resulting model is what we call fBm/D/1. Unfortunately, no simple, closed-form results have been found for this queuing model. However, Norros was able to prove the following approximating formula:

3.65 $$P_{tail}(h)\approx \exp \left(-\frac{1}{2am}{\left(\frac{(C-m)}{H}\right)}^{2H}{\left(\frac{h}{1-H}\right)}^{2-2H}\right)$$

where C is the speed at which the packets in the buffer are processed. In an infinite buffer (theoretical abstraction), P_tail(h) represents the prob-ability for the buffer length to exceed h bytes. (Regarded as a function of the variable h, the right-hand side is a Weibull distribution.)

There are two problems with this formula. From a theoretic stand point, the approximations made are nonsystematic (i.e., all in the same direction, either “≤” or “≥”), so Equation 3.65 is neither an upper bound nor a lower bound of the actual P_tail(h). The exponential could be either smaller or larger than P_tail(h) and there is no estimation of the tightness of the approximation. From an empirical standpoint, Equation 3.65 does not fit into any of our data series, most of the time deviating by orders of magnitude. Even assuming that for an exact fBm process Norros’ formula would be a good approximation, we still have to deal with the fact that our time series are not accurately described by fBm at low timescales.

3.3.3.3  M/G/∞ queue method

Let consider a queue M/G/∞ where the arrivals of the clients follow a Poisson distribution and the service time are obtained from a heavy-tailed distribution (with infinite dispersion). Then the process {X_t}, representing the number of clients in the system at the t moment, is asymptotically self-similar.

To exemplify this method, we present a mechanism for congestion avoidance in a packet-switched environment. If there are no resources reserved for a connection, the possibility of packets arriving at a node at a rate higher than the processing rate then appears.

Packets arriving on any of a number of input ports (A–Z) can be destined for the same output port. Although there could be, in general, many operations to be performed, they are broadly divided into two classes: serialization (i.e., the process of transmitting the bits out on the link) and everything else that happens inside the node (e.g., classification, rewriting of the TTL, and/or TOS fields in the IP header, routing table and access list look-ups, policing, shaping, etc.) called internal processing. Serialization delays are deterministic and determined by the speed of the output link. Internal processing delays are essentially random (a larger routing table requires a longer time to search) but advanced algorithms and data structures are doing, in general, a good job of maintaining overall deterministic performance. As a consequence, the only truly random part of the delay is the time a packet spends waiting for service.

The internal workings of a node are relatively well understood (by either analysis or direct measurement) and controllable. In contrast, the (incoming) traffic flows are neither completely understood nor easily controllable. The focus is therefore on accurately describing (modeling) these flows and predicting their impact on network performance.

The picture in Figure 3.3a is a simplified representation of a packet-switching node (router or switch), illustrating the mechanism of congestion for a simple M/D/1 queue model. The notation is easily explained: “M” stands for the (unique) input process, which is assumed Markovian, that is, a Poisson process; “D” stands for the service, assumed deterministic; and “1” means that only one facility is providing service (no parallel service).

It has been shown time and again that the high variability associated with LRD and SS processes can greatly deteriorate network performance, creating, for instance, queuing delays orders of magnitude higher than those predicted by traditional models. We make this point once again with the plot shown in Figure 3.3b: the average queue length created in a queuing simulation by the actual packet trace is widely divergent from the average predicted by the corresponding queuing model (M/D/1).

In view of our simple model of a node discussed earlier, hypotheses “D” and “1” seem reasonable, so we would like to identify “M” as the cause of the massive discrepancy. The question then arises: Is there an “fBm/D/1” queuing model? The answer is “almost,” which suggests that the self-similarity is mainly caused by user/application characteristics (i.e., Poisson arrivals of sessions, highly variable session sizes or durations) and is hence likely to remain a feature of network traffic for some time to come—assuming the way humans tend to organize information will not change drastically in the future.

Figure 3.3   A queue model (a) and the comparison of the real traffic with the simulated traffic (b).

3.3.3.4  Random midpoint displacement (RMD) method

The basic concept of the random midpoint displacement (RMD) algorithm is to extend the generated sequence recursively, by adding new values at the midpoints from the values at the endpoints. Figure 3.4 outlines how the RMD algorithm works.

Figure 3.5 illustrates the first three steps of the method, leading to generation of the sequence (d_3,1; d_3,2; d_3,3; and d_3,4). The reason for subdividing the interval between 0 and 1 is to construct the Gaussian increments of X. Adding offsets to midpoints makes the marginal distribution of the final result normal.

• Step 1. If the process X(t) is to be computed for time instance t between 0 and 1, then start out by setting X(0) = 0 and selecting X(1) as a pseudo-random number from a Gaussian distribution with mean 0 and variance $\text{Var}[X(1)]={\sigma }_0^2$
  . Then $\text{Var}[X(1)-1)(0)]={\sigma }_0^2$
  .
• Step 2. Next, X(1/2) is constructed as the average of X(0) and X(1), that is, X(1/2) = 1/2 [X(0) + X(1)] + d₁.
  

  Figure 3.4   Block scheme for the RMD method.

  

  Figure 3.5   The first three steps in the RMD method.

  The offset d₁ is a Gaussian random number (GRN), which should be multiplied by a scaling factor 1/2, with mean 0 and variance $S_1^2$
  of d₁. For $\text{Var}[X(t_2)-X(t_1)]=|t_2-t_1{|}^{2H}{\sigma }_0^2$
  to be true, for 0 ≤ t₁ ≤ t₂ ≤ 1, it must be required that $\text{Var}[X(1/2)-X(0)]+S_1^2={(1/2)}^{2H}{\sigma }_0^2$
  . Thus, $S_1^2={(1/2)}^{2H}(1-2^{2H-2}){\sigma }_0^2$
• Step 3. Reduce the scaling factor by $\sqrt{2}$
  , that is, now assume $\sqrt{8}$
  , and divide the two intervals from 0 and 1/2 and from 1/2 to 1 again. X(1/4) is set as the average 1/2 [X(0) + X(1/2)] plus an offset d_2,1, which is a GRN multiplied by the current scaling factor $\sqrt{1/8}$
  . The corresponding formula holds for X(3/4), that is, X(3/4) = 1/2 [X(1/2) + X(1)] + d_2,2 where d_2,2 is a random offset computed as before. So the variance $S_2^2$
  of d_2,* must be chosen such as follows:
  $$\text{Var}\left[X\left(\frac{1}{4}\right)-X(0)\right]=\left(\frac{1}{4}\right)\text{Var}\left[X\left(\frac{1}{2}\right)-X(0)\right]+S_2^2={\left(\frac{1}{2^2}\right)}^{2H}{\sigma }_0^2$$

  
  Thus, $S_2^2=(1/2^2)2^H(1-2^{2H-2}){\sigma }_0^2$
  .
• Step 4. The fourth step proceeds in the same manner: reduce the scaling factor by, that is, do $\sqrt{16}$
  . Then set
  $$\begin{array}{cc}\hfill X\left(\frac{1}{8}\right)=\frac{1}{2}\left(X(0)+X\left(\frac{1}{4}\right)\right)+d_{3,1};\hfill & \hfill X\left(\frac{3}{8}\right)=\left(\frac{1}{2}\right)\left(X\left(\frac{1}{4}\right)+X\left(\frac{1}{2}\right)\right)+d_{3,2}\hfill \end{array}$$

  
  $$\begin{array}{cc}\hfill X\left(\frac{5}{8}\right)=\frac{1}{2}\left(X\left(\frac{1}{2}\right)+X\left(\frac{3}{4}\right)\right)+d_{3,3};\hfill & \hfill X\left(\frac{7}{8}\right)=\left(\frac{1}{2}\right)\left(X\left(\frac{3}{4}\right)+X(1)\right)+d_{3,4}\hfill \end{array}$$

  
  In the formulas above, d_3,* is computed as a different GRN multiplied by the current scaling factor $\sqrt{16}$
  . The following step computes X(t) at t = 1/16, 3/16, …, 5/16 using a scaling factor again reduced by $\sqrt{2}$
  , and continues as indicated above. So the variance $S_3^2$
  of d_3,* must be chosen such that $\text{Var}[X(1/8)-X(0)]=(1/4)\text{Var}[X(1/5)-X(0)]+S_3^2=(1/2^3)2^H{\sigma }_0^2$
  . Thus, $S_3^2={(1/2^3)}^{2H}(1-2^{2H-2}){\sigma }_0^2$
  . The variance $S_n^2$
  of d_n*, therefore, yields $S_n^2={(1/2^n)}^{2H}(1-2^{2H-2}){\sigma }_0^2$
  . The main drawback of the method is the generation of a process, which is only approximate self-similar. In particular, the parameter H of the samples trajectories tend to overpass the target value for 0.5 < H < 0.75 and to be less than the target value for 0.75 < H < 1.0.

3.3.3.5  Wavelet transform-based method

This method (Abry and Veitch 1998) is based on the so-called multiresolution analysis (MRA), which consists of splitting the sequence into a (low pass) approximation and a series of details (high pass). First, the coefficients corresponding to the wavelet transform of fBm. Then, these are used in an inverse wavelet transform in order to obtain fBm samples. The disadvantage is that it is only approximate, because generally the wavelet coefficients are not correlated. Nevertheless, the problem can be compensated in the case of wavelets with a sufficient number of moments reduced to zero (i.e., the degree of the polynomial for that the internal scalar product with the wavelet is zero).

For some time after the discovery of scaling in traffic, it was debated as to whether the data was indeed consistent with self-similar scaling, or that the finding was merely an artifact of poor estimators in the face of data polluted with nonstationarities. The introduction of wavelet-based estimation techniques to traffic helped greatly to resolve this question, as they convert temporal LRD to SRD in the domain of the wavelet coefficients, and simultaneously eliminate or reduce certain kinds of nonstationarities. It is now widely accepted that scaling in traffic is real, and wavelet methods have become the method of choice for detailed traffic analysis. Another key advantage of the wavelet approach is its computational complexity (memory complexity is also in an online implementation), which is invaluable for analyzing the enormous data sets of network-related measurements. However, even wavelet methods have their problems when applied to certain real-life or simulated traffic traces. An important ruleof-thumb continues to be to use as many different methods as possible for checking and validating whether or not the data at hand is consistent with any sort of hypothesized scaling behavior, including self-similarity.

Wavelet analysis is a joint timescale analysis. It replaces X(t) by a set of coefficients, such as d_X(j, k), j, k ∈ Z where 2^j denotes the scale, and the 2^jk instant about which the analysis is performed. In the wavelet domain, Equation 3.57 is replaced by Var[d_X(j, k)] = c_fC2^(1−β)j, where the role of m is played by scale, of which j is the logarithm, where c_f is the frequency domain analog of c_γ, and C is independent of j. The analog of the variance– time diagram, which is the estimates of log(Var[d_X(j, ⋅)]) against j, called the logscale diagram. It constitutes an estimate of the spectrum of the process in log–log coordinates, where the low frequencies correspond to large scales, appearing in the right in the figures below. The global behavior of data, as a function of scale, can be efficiently examined in the logscale diagram, and the exponent estimated by a weighted linear regression over the scales where the graph follows a straight line. What constitutes a straight line can be judged using the confidence intervals at each scale, which can be calculated or estimated. The important fact is that the estimation is heavily weighted toward the smaller scales, where there is much more data.

3.3.3.6  Successive random addition (SRA)

Another alternative method for the direct generation of fBm process is based on the successive random addition (SRA) algorithm. The SRA method uses the midpoints like RMD, but adds a displacement of a suitable variance to all of the points to increase stability of the generated sequence.

Figure 3.6   Block scheme of the SRA method.

Figure 3.6 shows how the SRA method generates an approximate self-similar sequence. The reason for interpolating midpoints is to construct Gaussian increments of X, which are correlated. Adding offsets to all points should make the resulted sequence self-similar and of normal distribution.

The SRA method consists of the following four steps:

• Step 1. If the process {X_t} is to be computed for time instance t between 0 and 1, then start out by setting X₀ = 0 and selecting X₁ as a pseudorandom number from a Gaussian distribution with mean 0 and variance $\text{Var}[X_1]={\sigma }_0^2$
  . Then $\text{Var}[X_1-X_0]={\sigma }_0^2$
  .
• Step 2. Next, X_1/2 is constructed by the interpolation of the midpoint, that is, X_1/2 = (1/2) (X₀ + X₁).
• Step 3. Add a displacement of a suitable variance to all of the points, that is, X₀ = X₀ + d_1,1, X_1/2 = X_1/2 + d_1,2, and X₁ = X₁ + d_1,3. The offsets d_1,* are governed by fractional Gaussian noise. For $\text{Var}[X(t_2)-X(t_1)]=|t_2-t_1{|}^{2H}{\sigma }_0^2$
  to be true, for any t₁, t₂, 0 ≤ t₁ ≤ t₂ ≤ 1, it must be required that Var[X / −X ] = (1/4) $\text{Var}[X_1-X_0]+2S_1^2={(1/2)}^{2H}{\sigma }_0^2$
  . Thus, $S_1^2=(1/2){(1/2^1)}^{2H}(1-2^{2H-2}){\sigma }_0^2$
  .
• Step 4. Next, Steps 2 and 3 are repeated. Therefore $S_n^2=(1/2){(1/2^n)}^{2H}(1-2^{2H-2}){\sigma }_0^2$
  , where ${\sigma }_0^2$
  is an initial variance and 0 < H < 1. Using the above steps, the SRA method generates an approximate self-similar fBm process.

3.5.1.1  Network traffic

Network traffic is obviously driven by applications. The most common applications that come to mind might include Web, e-mail, peer-to-peer file sharing, and multimedia streaming. A traffic study of Internet backbone traffic showed that Web traffic is the single largest type of traffic, on the average more than 40% of total traffic (Fraleigh et al. 2003). Although on a minority of links, peer-to-peer traffic was the heaviest type of traffic, indicating a growing trend. Streaming traffic (e.g., video and audio) constitute a smaller but stable portion of the overall traffic. E-mail is an even smaller portion of the traffic. Therefore, traffic analysts have attempted to search for models tailored to the most common Internet applications (Web, peerto-peer, and video).

3.5.1.2  Web traffic

The World Wide Web is a familiar client–server application. Most studies have indicated that Web traffic is the single largest portion of overall Internet traffic, which is not surprising considering that the Web browser has become the preferred user-friendly interface for e-mail, file transfers, remote data processing, commercial transactions, instant messages, multimedia streaming, and other applications.

An early study focused on the measurement of self-similarity in Web traffic (Crovella and Bestavros 1997). Self-similarity had already been found in Ethernet traffic. Based on days of traffic traces from hundreds of modified Mosaic Web browser clients, it was found that Web traffic exhibited self-similarity with the estimated Hurst parameter H in the range 0.7–0.8. More interestingly, it was theorized that Web clients could be characterized as on/off sources with heavy-tailed on periods. In previous studies, it had been demonstrated that self-similar traffic can arise from the aggregation of many on/off flows that have heavy-tailed on or off periods. Heavy tails implies that very large values have nonnegligible probability. For Web clients, on periods represent active transmissions of Web files, which were found to be heavy-tailed. The off periods were either “active off,” when a user is inspecting the Web browser but not transmitting data, or “inactive off,” when the user is not using the Web browser. Active off periods could be fit with a Weibull probability distribution, while inactive off periods could be fit by a Pareto probability distribution. The Pareto distribution is heavy-tailed. However, the self-similarity in Web traffic was attributed mainly to the heavy-tailed on periods (or essentially, the heavy-tailed distribution of Web files). For example, let consider that a Web page consists of a main object, which is an HTML document, and multiple in-line objects that are linked to the main object (see Figure 3.8). Both main object sizes and in-line object sizes are fit with (different) lognormal probability distributions. The number of in-line objects associated with the same main object follow a gamma probability distribution. In the same Web page, the intervals between the start times of consecutive in-line objects are fit with a gamma probability distribution.

Figure 3.8   ON/OFF model with Web page.

TCP flow models seem to be more natural for Web traffic than “page based” consisting of a main object linked to multiple in-line objects.

In this example, the requests for TCP connections originate from a set of Web clients and go to one of multiple Web servers. A model parameter is the random time between new TCP connection requests. A TCP connection is held between a client and server for the duration of a random number of request–response transactions. Additional model parameters include sizes of request messages, server response delay sizes of response messages, and delay until the next request message. Statistics for these model parameters are collected from days of traffic traces on two high-speed Internet access links.

3.5.1.3  Peer-to-peer traffic

Peer-to-peer (P2P) traffic is often compared with Web traffic because P2P applications work in the opposite way from client–server applications, at least in theory. In actuality, P2P applications may use servers. P2P applications consist of a signaling scheme to query and search for file locations, and a data transfer scheme to exchange data files. While the data transfer is typically between peer and peer, sometimes the signaling is done in a client–server way. Since P2P applications have emerged only in the last decade, there are a few traffic studies to date. They have understandably focused more on collecting statistics than the harder problem of stochastic modeling. An early study examined the statistics of traffic on a university campus network (Saroiu et al. 2002).

The observed P2P applications were Kazaa and Gnutella. The bulk of P2P traffic was AVI and MPEG video, and a small portion was MP3 audio. With this mix of traffic, the average P2P file was about 4 MB, three orders of magnitude larger than the average Web document. A significant portion of P2P was dominated by a few very large objects (around 700 MB) accessed 10–20 times. In other words, a small number of P2P hosts are responsible for consuming a large portion of the total network bandwidth.

3.5.1.4  Video

The literature on video traffic modeling is vast because the problem is very challenging. The characteristics of compressed video can vary greatly depending on the type of video that determines the amount of motion (from low-motion videoconferencing to high motion broadcast video) and frequency of scene changes; and the specific intraframe and interframe video coding algorithms. Generally, interframe encoding takes the differences between consecutive frames, compensating for the estimated motion of objects in the scene. Hence, more motion means that the video traffic rate is more dynamic, and more frequent scene changes is manifested by more “spikes” in the traffic rate. While many video coding algorithms exist, the greatest interest has been in the Moving Picture Coding Experts Group (MPEG) standards. Video source rates are typically measured frame by frame (which are usually 1/30 s apart). That is, X_n represents the bit rate of the nth video frame. Video source models address at least these statistical characteristics: the probability distribution for the amount of data per frame; the correlation structure between frames; times between scene changes; and the probability distribution for the length of video files. Due to the wide range of video characteristics and available video coding techniques, it has seemed impossible to find a single universal video traffic. The best model varies from video to video file. A summary of research findings was presented by Chen (2007) and includes the following:

• The amount of data per frame has been fit to log-normal, gamma, and Pareto probability distributions or their combinations.
• The autocorrelation function for low-motion video is roughly exponential, while for broadcast video is more complicated, exhibiting both SRD at short lags and LRD at much longer lags.
• Time series, Markov modulated, and wavelet models are popular to capture the autocorrelation structure.
• The lengths of scenes usually follow Pareto, Weibull, or Gamma probability distributions.
• The amount of data in scene change frames are uncorrelated and follow Weibull, Gamma, or unknown probability distributions.
• The lengths of streaming video found on the Web appear to have a long-tailed Pareto probability distribution.

3.5.2.1  TCP flows with congestion avoidance

TCP senders follow a self-limiting congestion avoidance algorithm that adjusts the size of a congestion window according to the inferred level of congestion in the network. The congestion window is adjusted by so-called additive increase multiplicative decrease (AIMD). In the absence of congestion when packets are acknowledged before their retransmission time-out, the congestion window is allowed to expand linearly. If a retransmission timer expires, TCP assumes that the packet was lost and infers that the network is entering congestion. The congestion window is collapsed by a factor of a half. Thus, TCP flows generally have the classic “saw tooth” shape. Accurate models of TCP dynamics are important because of TCP’s dominant effect on Internet performance and stability. Hence, many studies have attempted to model TCP dynamics (a significant synthesis is presented in Altman et al. 2005). Simple equations for the average throughput exist, but the problem of characterizing the time-varying flow dynamics is generally difficult due to interdependence on many parameters such as the TCP connection path length, round-trip time, and probability of packet loss, link bandwidths, buffer sizes, packet sizes, and number of interacting TCP connections.

3.5.2.2  TCP flows with active queue management

The TCP congestion avoidance algorithm has an unfortunate behavior where TCP flows tend to become synchronized. When a buffer overflows, it takes some time for TCP sources to detect a packet loss. In the meantime, the buffer continues to overflow, and packets will be discarded from multiple TCP flows. Thus, many TCP sources will detect packet loss and slow down at the same time. The synchronization phenomenon causes underutilization and large queues.

Active queue management schemes, mainly random early detection (RED) and its many variants, eliminate the synchronization phenomenon and thereby try to achieve better utilization and smaller queue lengths. Instead of waiting for the buffer to overflow, RED drops a packet randomly with the goal of losing packets from multiple TCP flows at different times. Thus, TCP flows are unsynchronized, and their aggregate rate is smoother than before. From a queuing theory viewpoint, smooth traffic achieves the best utilization and shortest queue.

In conclusion, the impacts of the TCP flow control mechanisms on the traffic are very strong. TCP-type feedback appears to have the effect of modifying the self-similarity traffic behavior, but it neither generates it nor eliminates it.