Generation of synthetic video traffic using time series

Simulation Modelling Practice and Theory 75 (2017) 127–145

Contents lists available at ScienceDirect

Simulation Modelling Practice and Theory journal homepage: www.elsevier.com/locate/simpat

Generation of synthetic video traffic using time series Christos Katris∗, Sophia Daskalaki Department of Electrical & Computer Engineering, University of Patras, Rio 26504 Greece

a r t i c l e

i n f o

Article history: Received 9 June 2016 Revised 5 February 2017 Accepted 13 April 2017

Keywords: Synthetic video traffic Time series models FARIMA FARIMA-GARCH Non-normality FARIMA-to-LN

a b s t r a c t Generating traffic has always been an important part of network simulations but has turned to an even more challenging task with modern networks. The statistical properties of the input stochastic processes traced in the networks used all along Information Era turned out to be complicated and difficult to reproduce. Taking into account successful efforts in modeling Internet traffic with FARIMA time series models, this paper attempts to extend their applicability and employ them to generate synthetic video traffic. It is known that FARIMA can model both the Short Range (SRD) and Long Range Dependence (LRD) existing in video traffic; however the traces it produces fail to describe correctly the moments (mean, standard deviation, skewness, kurtosis) of the distribution behind the data. Since an efficient traffic generator should capture both the statistical properties and queuing behavior of video traffic we experiment with models such as FARIMA with Student’s t errors and FARIMA-GARCH with Normal and Student’s t errors, improving somewhat the accuracy of the generated traffic. Furthermore, the paper suggests the projection of the traces generated by a FARIMA model to values of a Lognormal distribution. It is shown that such a methodology produces synthetic traces that can emulate very closely the behavior of real traces. In order to quantify closeness the generated traces are fed into a simple FIFO queuing system with finite buffers, where loss probability is calculated and compared to that experienced by the corresponding real traces. Using five different real traces, MPEG-4 or H.263, it is shown that the proposed methodology produces traffic generators that can capture satisfactorily several statistical properties of the real traffic and also its queuing behavior for a wide range of buffer sizes and service rates. © 2017 Elsevier B.V. All rights reserved.

1. Introduction Generation of synthetic traffic is an important procedure quite necessary for performance evaluation when designing a network and for network management. Its necessity has become even more pressing due to the remarkable growth of video applications and multimedia users reported during the recent years. When reliable, synthetic traffic generators can be used with simulations to resolve problems in multimedia traffic management, such as bandwidth allocation, design of protocols and admission controls [1], while testing real networks for their performance can be expensive and sometimes prohibitive. For this reason it is preferable to use mathematical models and/or simulations to evaluate the performance of a network. In general, statistical models are considered better choice than trace-driven simulations since they can shed light to traffic characteristics, they are stochastic and not static, and can describe different situations when tuning the parameters in the models. ∗

Corresponding author. E-mail addresses: [email protected] (C. Katris), [email protected] (S. Daskalaki).

http://dx.doi.org/10.1016/j.simpat.2017.04.005 1569-190X/© 2017 Elsevier B.V. All rights reserved.

128

C. Katris, S. Daskalaki / Simulation Modelling Practice and Theory 75 (2017) 127–145

In order to generate synthetic video traffic using time series one should consider models that are suitable for this kind of traffic, such as FARIMA (Fractional ARIMA), since these models capture properties like SRD (Short Range Dependency) and LRD (Long Range Dependency) that normally video traffic appear to carry. According to Alheraish [2] a successful video traffic model should be able to capture the statistical characteristics of a real trace (mean, standard deviation, coefficient of variation, mode and autocorrelation) and also its dynamical behavior, e.g. to expose similar queuing behavior. It should be as simple as possible, generate traffic with low computational complexity and if possible characterize a wide range of video sources. According to Casilari et al. [3] approximating the probability density function (pdf) of the underlined random process is of utmost importance for matching the queuing behavior of a real traffic trace, while SRD is of secondary importance and LRD of even less. The method proposed in this paper satisfies relatively well the above criteria and certainly improves proximity to real traffic compared to the classical FARIMA models. The main scope of the paper is to present a new procedure that is based on time series models and can be used for generating synthetic traffic for full motion video. Several VBR source models that have been used in the past are based on Autoregressive (AR) processes and a comprehensive description can be found in [4], along with a thorough and most recent survey on models for video traffic. Here we briefly refer only to those that have been used for video traffic generation. In general the AR approach is simple and requires few parameters. For full motion video, however, which could comprise frames of different types (I, B, or P), single AR models are not suitable and usually a combination of AR processes with different residual distributions for each frame type is required. As a consequence frame-based [5,6], scene-based and nested AR models [7] have been adopted. Of course this increases complexity to modeling and formulation. Other alternative variations of AR models include the DAR (Discrete Autoregressive) model of order p [8,9], GAR (AR model with Gamma residual process), GBAR (Gamma-Beta Autoregressive model) [10] and the General AR model that also includes GACS (Gaussian Autoregressive and Chi-Square) [11,12]. The projected AR (PAR) models form yet another interesting category of models for generating video traffic. They attempt to preserve autocorrelation and at the same time to fit the frame-size histogram. Previous work on such traffic models uses different types of AR processes according to the type of video and encoding scheme [13]. They appear to capture in general the autocorrelation behavior of a video and estimation of their coefficients from empirical data is relatively simple. On the other hand, it is not possible to find a single AR model that can capture different statistical characteristics and also there is no single video model that is suitable for all video sequences and for all purposes. Moreover, self-similar models have been used successfully for video traffic modeling and generation. FARIMA models are the most widely used, individually [14,15] and in more composite frameworks [16,17,3]. Casilari et al. [3] considered briefly projections of the FGN (Fractional Gaussian Noise) and FARIMA models to a desired distribution. Apart from the time series models there has been a long standing effort for modeling traffic using Markov and Markov modulated models, most of them quite successfully for some type of video [e.g. [18,19]]. Their level of complexity is considered much higher than the level required from the AR family of models. Lastly, wavelets have formed another group of successful models for video traffic. Their success has been reported in a number of articles [e.g. [20,21]]. The ultimate goal of our approach is to provide a general framework, which can be adapted to several circumstances and is able to generate traffic that mimics real video traces, using only a certain time unit and very little information of the empirical traffic, such as the mean and standard deviation. In the core of our approach is the dependence structure (long and short range) that exists in video traffic, and for this reason we employ FARIMA models. On the other hand, for capturing the queuing behavior of real traffic, the marginal distribution of the underlying process should be matched. For this reason, we consider a 2-stage approach where an initially fitted FARIMA model to a video trace is transformed to a data series with a pre-determined target distribution. It builds upon the PAR types of models using FARIMA modeling in conjunction with the Lognormal distribution, which requires only two parameters and its CDF is easily invertible. With this approach we can accurately capture the marginal distribution of the traffic. Although this approach is superior to other FARIMA and ARIMA alternatives in capturing queuing performance of real traffic, we also have to assess its ability to preserve the SRD and LRD property of the FARIMA model. We perform this assessment experimentally. The classical FARIMA model is usually very successful in capturing the LRD and SRD properties of a real traffic trace. However our experimentation showed that it cannot capture its marginal distribution, when it is far from normal, neither its non-linear dependencies, when they are present. To overcome this we considered as traffic generators FARIMA models with Student’s t and Normal Inverse Gaussian error distributions, and also FARIMA-GARCH models with Normal, Normal Inverse Gaussian and Student’s t errors. A number of these models achieved slight improvements for some traces, but none achieved a clear and definite improvement. As an alternative we develop a two-step traffic generating procedure which combines FARIMA models with Lognormal distribution through a transformation of percentiles. Such a procedure first embodies the desired SRD and LRD characteristics with the help of a FARIMA model and then attempts to approximate the pdf of a recommended distribution, the Lognormal in this case. The method is analogous to ARTA (AR To Anything) process [22,23] which combines an AR model with a desired marginal distribution and generates a series that displays the autocorrelation of the AR model and at the same time achieves a desired marginal distribution. The method is tested with simulated traces which include different characteristics, and then applied to real traces successfully compared with the aforementioned models. Next, we describe analytically the procedure, argue on its efficacy and set up a queuing framework for evaluating its performance. The evaluation is performed experimentally and shows that transforming a FARIMA generated trace to variates of a Lognormal distribution matches sufficiently well real video traces. Section 2 presents the time series models that will

C. Katris, S. Daskalaki / Simulation Modelling Practice and Theory 75 (2017) 127–145

129

be used in this paper, i.e. FARIMA and FARIMA-GARCH. Section 3 gives the details for the proposed procedure which can be used either to generate video traffic when significant amount of information is known about a video trace but also when limited information is given. In Section 4 a performance evaluation scheme based on a FIFO single server queuing system is proposed for the comparison of the produced traces. Section 5 renders an extensive experimentation of the procedure using five different video traces and lastly Section 6 summarizes and concludes this work. 2. Synthetic traffic with time series models In this section we present the time series models that will be used for generating traffic. It is assumed that traffic is expressed with a random process {Xt , t = 1, 2, . . .} which represents the traffic load that arrives at time unit t. We consider discrete-time time series models and focus on those that match the statistical characteristics of video traffic, i.e. FARIMA and FARIMA-GARCH. 2.1. FARIMA models The Fractional ARIMA (FARIMA) models introduced by Granger and coworkers [24,25], are extensions of ARIMA and have the ability to describe both short and long range dependencies of the data. FARIMA models have been used previously to describe internet traffic in general and video traffic in particular. In [26] FARIMA models were used to model broadband traffic and in several situations they were found to improve pure LRD models. In [27] FARIMA models were found to be suitable for estimating the queuing performance under heavy-load conditions. In [28] a FARIMA based methodology for bandwidth allocation was reported to be more successful compared to other models. Similarly, in [29] FARIMA models were used as basis in a bandwidth allocation scheme for different types of traffic, including JPEG and MPEG traffic. Lastly, in [1] and [3] FARIMA models are compared with White Noise, Autoregressive (AR) and Fractal Gaussian Noise (FGN), and found to be superior in modeling VBR traffic. The formulation of a FARIMA time series model is given as:

p (L )(1 − L )d (Xt ) = q (L )εt , where L is the lag operator, p (L ) = 1 − ϕ1 L − . . . − ϕ p L p , q (L ) = 1 + θ1 L + . . . + θq Lq , (1 − L )d =

  d j

(−1 ) =

∞

j=0

  d j

(1)

(−1 ) j L j and

(−d+ j ) (−d )( j+1 ) .

Long memory is displayed by the process, when 0 < d < 0.5 . The classical FARIMA model assumes that the error terms ɛt are normally distributed with zero mean and constant variance σ 2 . However, since studies with Internet traffic datasets suggest mostly non-symmetric leptokurtic distributions, for several applications the classical model carries an obvious weakness. For this reason, we consider FARIMA models with Student’s t innovations, i.e. the error terms are assumed to follow Student’s t distribution with zero mean, constant variance σ 2 and v degrees of freedom, where v > 2. The pdf of the location-scale version of Student’s t distribution is known to be:



 v+1 



( x − α )2 f (x; a, β , v ) =  v 1 + βv βvπ  2 2

−( v+1 2 )

where α , β and v are the location, scale and shape parameters, respectively, and  (.) is the Gamma function. The mean of the distribution is the location parameter α and the variance is given as (vβv . We should note here that for this study we −2 ) additionally experimented with Normal Inverse Gaussian distribution for the error terms but experienced large deviations for several traffic traces, so there is no need reporting their performance. The fitting procedure for a FARIMA model to a dataset is not trivial in general and the one followed here is described in [30]. The estimation of the parameters (except d) is realized through Maximum Likelihood (MLE), however a crucial step is to decide upon the value of d, the fractional parameter, and for this paper it is estimated using the Geweke and PorterHudak estimator [31]. When the estimator fails to achieve a value between 0 and 0.5, then we conventionally use the value 0.25. 2.2. FARIMA-GARCH models The FARIMA-GARCH models are extensions of the FARIMA models, in the sense that they allow heteroskedasticity for the errors and in some cases may achieve better results than plain FARIMA. These models assume conditionally heteroskedastic errors, quite a rational assumption for video traffic, since in a video there are scenes that create lumps of heavy traffic (such as action scenes, explosions etc.) and others where traffic is much lighter. FARIMA-GARCH models have been used for the prediction of inflation in [32] and when variance cannot be assumed as stationary, we expect that their usage for generating synthetic traffic will improve both modeling and queuing performance. Examples of articles where linear models such as ARIMA are combined with the non-linear GARCH models for internet traffic modeling are [33], where the author concludes that the seasonal AR-GARCH models are superior over the plain seasonal ARIMA models in forecasting accuracy;

130

C. Katris, S. Daskalaki / Simulation Modelling Practice and Theory 75 (2017) 127–145

[34], where ARIMA-GARCH is found to have larger forecasting accuracy than FARIMA; and [35] where it is stated that in that study ARIMA-GARCH model was superior than Seasonal ARIMA (SARIMA) model for forecasting HTTP traffic. In this paper, we consider FARIMA-GARCH models, because they have the ability to capture SRD and LRD and because of conditional heteroskedasticity. A FARIMA-GARCH model follows the equation:

 p (L )(1 − L )d (Xt − μ ) = q (L )εt where

εt | t−1 ∼ probabilit y dist ribut ion ( p1,..., pk ) and σt2 = ω +

(2)

q

2 αi εt−1 + i=1

p i=1

2 βi σt−1 .

Comparing with the FARIMA model in Eq. (1), the new element here is the second line, where it implies that the error term has a conditional PDF with parameters ( p1,..., pk ). Most commonly it is assumed that the distribution conditional on the information set t−1 is normal with zero mean and σt2 variance. In the last equation, the conditional variance follows the GARCH model of [36]. The fitting procedure we use in this paper is the following. First, we consider the same order and same fractional parameter d as for the corresponding FARIMA process. Additionally, for the conditional variance we consider a GARCH(1,1) model, which is generally accepted to be a suitable approximation for most applications. These two steps are performed for errors that follow Normal or Student’s t distributions. After the fitting procedures are finalized the resulting models are used for traffic generation through simulation. For the first part of our experimental work, shown in Section 5, we attempt generating traffic using FARIMA and FARIMAGARCH models with Normal and Student’s t errors. The results indicate that FARIMA with t errors prevail and the generated traces follow relatively closer the behavior of the real traces. As a next step and in order to approximate even further we introduce a 2-step methodological procedure that in addition involves the marginal distribution of Xt . 3. Proposed approach for generating video traffic The 2-step approach to be introduced here aims at generating streams of synthetic traffic that behave similarly to some existing video stream and for this we attempt to approximate satisfactorily both the dependence structure observed in the data and its marginal PDF. First, we create a sequence of values with the help of a time series model which is suitable for describing the dependence structure found in an actual video stream; and then through a statistical transformation, which is equivalent to the Inverse Sampling method, we project the simulated values to a suitable distribution so to approximate the marginal distribution hypothesized for the actual data. The rationale behind this procedure stems from the fact that the initial simulated data from the time series model will somehow exhibit the dependence structure traced in the actual video stream; however, as it became clear from our study, they will not normally fit well the PDF of the real stream. In order to handle such a deficiency, the simulated data are projected to a new corresponding trace carrying the desired distribution. The projection is realized using the percentile ranks of the data simulated previously and then through the inverse CDF that corresponds to some distribution assumed for the real traffic trace. As it turns out this procedure is capable of creating streams that exhibit approximately a desired dependence structure and at the same time fit well to some desired distribution so it can finally behave very closely to the original trace. 3.1. Description of the procedure Given a real video traffic trace let our goal be the generation of a synthetic trace of length N that behaves similarly to a real trace. Step 1. Simulate N observations from a suitable time series model. We consider a time series model that can describe successfully the dependence structure of the data. Such a model could be either a FARIMA or FARIMA-GARCH and initially it is applied as generator for the synthetic traffic (yi , i = 1, . . . , N). Based on theory it is expected that the simulated data will exhibit a dependence structure similar to that of the real observations (i.e. SRD and LRD dependencies). Regarding the distribution however, when we use the classical FARIMA model the generated data will have approximately a normal distribution: (yi ∼ N(μ, σ 2 )). Step 2. Transform simulated observations from step 1 to a target distribution. Given a “target” marginal distribution D, the second step of the traffic generation is to transform the simulated data to follow its PDF and at the same time retain the original dependence structure. The proposed transformation is realized through percentiles. More specifically, the probabilities pi  s that correspond to theyi  s (traffic from step 1) are all calculated using the Normal distribution ( pi = FN (yi )). Then using the inverse cumulative function of D the percentiles xi = FD−1 ( pi ) that correspond to the pi  s are calculated. The values xi produced form a new trace of synthetic traffic distributed according to D with mean μ ˆ D and standard deviation σˆ D and carry approximately the same dependence structure as the yi  s, both essential ingredients for capturing the dynamics of the behavior of the real traces in a queuing situation.

C. Katris, S. Daskalaki / Simulation Modelling Practice and Theory 75 (2017) 127–145

131

3.2. Selecting the target distribution Based on previous studies the marginal distribution of internet traffic cannot be described by a normal distribution. For video traffic, and in general for network traffic, the lognormal has been suggested instead [4,37]. In order to evaluate the suitability of Lognormal for our case, the empirical distribution is constructed first using kernels [38,39]. The kernel  x−x density estimate for the empirical function of a set of data is given by fˆh (x ) = 1n ni=1 1h K ( h i ), where K(.) is the kernel, a non-negative function that has mean zero and integrates to 1, and h is a smoothing parameter called bandwidth. There are several kernel functions such as uniform, triangular, Epanenchnikov and normal. For this paper we used the normal (Gaussian) kernel and the univariate plug-in selection method of [40] for the selection of the bandwidth. In Section 5.3 the empirical distributions of five different real video traces will be calculated and compared to Lognormal models to conclude that among different choices the Lognormal distribution prevails; this is due to its goodness of fit but also because it requires only two parameters that are easily estimated with the knowledge of a targeted mean and variance. 3.3. Evaluation of dependence structure The 2-step methodology previously described is based on the assumption that the traces generated from such a procedure will retain, at some significant level, the dependence structure originally held from the time series. The dependence structure here refers mainly to short and long range dependencies, i.e. autocorrelation and Hurst exponent, and to possible non-linearities. In order to test such a premise, in Section 5 we proceed with an additional experimentation. At first, data are simulated from ARIMA, FARIMA, ARIMA-GARCH and FARIMA-GARCH models. It is expected that the simulated data carry the SRD property and this can be diagnosed with a Ljung-Box test indicating significant autocorrelation. The Hurst exponent (H), on the other hand, can indicate the presence of LRD. It is a measure of long memory and takes values from 0 to 1. If H ∈ (0, 0.5) the series is expected to be antipersistent (i.e. high values tend to follow low values and vice versa), if H ∈ (0.5, 1) the series is expected to exhibit long-memory (i.e. high values tend to follow high values for some time into the future) and if H = 0.5 then absence of long memory is expected. There are a few different methods to estimate the H exponent in a time series, however there is no general agreement about the most accurate and/or suitable method. In this paper we will employ the Aggregated Variance method [41]. However, towards the end of our experimental work in Section 5, we additionally consider other alternatives for the calculation of Hurst exponent, including the more recent and according to studies accurate wavelet method [42]. Lastly, for the property of non-linear dependencies one can use the White Neural Network test. With these three tools it is possible to trace the dependencies existing in each one of the simulated series. Then using the inverse cumulative function of a Lognormal model the aforementioned transformation can be applied to the simulated data and thus generate traces that follow Lognormal distributions of desired parameters. Since it is not certain at which level the dependence structure is carried over, the transformed data will be re-evaluated using the same tools as previously. 3.4. Application to a new video trace The procedure presented in Section 3.1 assumes availability of a real trace which is actually utilized to fit the initiating time series model. In case of an unknown video trace this step will not be possible; however the proposed procedure can still be used as long as there is an estimate or a guess for its mean and variance. Given there is no trace to fit the model, a FARIMA(1,1) can be considered as reference model, with coefficients p1 = 0.5, q1 = 0.25 and fractional parameter d = 0.3. This value of d corresponds to a Hurst exponent of approximately H = 0.8 and establishes long range dependencies in the process, quite a realistic assumption for real video traces. At the second step of the procedure an additional problem is related to the parameter estimation (μ ˆ logyi , σˆ logyi ). Since it is assumed that there is an estimate (or at least a target value) for the mean and standard deviation of the trace, then it is possible to estimate the Lognormal parameters through the Bayesian Monte Carlo Markov Chain (MCMC) approach. In Section 5.5 the modified procedure just described is applied to generate traffic using only the means and standard deviations of several traces. The procedure’s efficacy is then evaluated by comparing the performance of the simulated and real traces in the queuing framework described in the next section. 4. Empirical performance evaluation of the approach In traffic generation it is well understood that similarity of traces can be achieved when their corresponding descriptive statistics, PDF and time-dependent characteristics are matched, since they all influence somehow the performance of a trace. However, since it is not easy to set the level at which all such characteristics should be matched throughout the empirical analysis we will follow a 3-fold comparison. One refers to the descriptive statistics, i.e. mean, standard deviation, skewness, kurtosis, coefficient of variation (CV), the other to certain time-dependent characteristics, i.e. stationarity, autocorrelation, LRD, non-linearity, and the third to the dynamical behavior of a trace which is formed by the video sequences. While the first two are pretty standard to handle, for the third we propose a queuing framework where the performance of traffic traces can be evaluated and compared. If the queuing performance of the synthetic traffic in this queue is close to the performance of the traffic it imitates the generating procedure is validated to be successful.

132

C. Katris, S. Daskalaki / Simulation Modelling Practice and Theory 75 (2017) 127–145 Table 1 Overview of datasets.

Trace Tr1 Tr2 Tr3 Tr4 Tr5

-

Silence of the Lambs Mr. Bean Star Trek Die Hard III The Firm

Type

Length (in 12-frame aggregates)

VBR / H.263 VBR / H.263 VBR / H.263 MPEG4 (HQ) VBR / H.263

40 0 0 40 0 0 40 0 0 40 0 0 40 0 0

Table 2 Descriptive statistics of real traces. Trace

Mean

Sd

Skewness

Kurtosis

CV

Tr1 Tr2 Tr3 Tr4 Tr5

32.108 31.608 17.640 42.528 14.309

31.343 16.977 12.799 22.150 10.915

3.256 1.922 2.518 1.363 2.143

19.438 10.453 14.158 5.325 10.112

0.976 0.537 0.726 0.521 0.763

Let us consider a single server FIFO queue with deterministic service rate and finite buffer. The queue is fed with traffic which is either a real trace or traffic generated from one of the procedures described previously. At each time period some traffic load (in kBytes) arrives and requests bandwidth. If the load exceeds availability (service capacity), the excess is placed in the buffer; however if the buffer size is also exceeded then loss will be realized. The traffic that should be served at time t (TRt ), at a FIFO basis, is the buffered load up until time t-1 (QTt -1 ) plus the new traffic that arrives during time period t (NTt ).

T Rt = N Tt + Q Tt−1 The performance measure that will be used for the queuing system under study is the loss probability: Pr[T Rt > SC + BS], where SC and BS are the service capacity and buffer size, respectively. The evaluation will be performed in two ways. Firstly the generators will be compared for their ability to capture the queuing behavior of real traffic when service rate is fixed and buffer sizes vary, and secondly when the buffer size is set to some fixed value and the service rates vary. Using these two schemes for comparison the traces generated by the different models will be compared for their performance with the actual trace and among themselves. For each traffic model the goal is to generate traces that are as close as possible to the real traffic traces, where closeness refers not only to the descriptive statistics, but also to the dependence structure characteristics, i.e. SRD and LRD, and also its queuing behavior for a wide range of service rates and buffer sizes. 5. Empirical analysis with real video traces In this section, the methodology described previously is realized to generate traffic that successfully emulates the behavior of real video traces. The focus is on capturing, as much as possible, the key statistical characteristics found in a collection of real traces, so that in a typical queuing situation the generated synthetic traces would perform similarly to their corresponding real traces. 5.1. Description of the real traces The data which are used throughout this analysis consist of five traces, four VBR with H.263 coding and one MPEG-4 of high quality. The five datasets are excerpts from publicly available video traces in TU-Berlin [43], while the aggregation level was decided to be the same for all of them and equal to kBytes per 12 frames (approximately a GOP level). Tables 1 and 2 display an overview of the datasets and their descriptive statistics, respectively. In Table 1 we give the names of the videos where the traces came from, the types and coding schemes, and the length considered for each dataset. In Table 2, for each trace we report the mean, standard deviation, skewness, kurtosis and coefficient of variation (CV). All datasets indicate clearly right skewed and leptokurtic distributions, with Tr4 to display the smallest discrepancy from symmetry and possibly normality. It is also interesting to note that the CV, which is a measure of relative dispersion, varies among the traces and it is possible that larger values indicate increased difficulty in modeling the trace. Moreover, Table 3 presents certain time-dependent properties observed for the five traces, i.e. stationarity, autocorrelation, long range dependence and non-linearity. According to the ADF (Augmented Dickey-Fuller) and Ljung-Box tests all traces can be assumed stationary and also display strong autocorrelation. Additionally all traces indicate the existence of LR dependency, with the estimated Hurst exponents to be in the range (0.5, 1.0) or, as in the case of Tr4, barely over 1. Finally, using the White NN test [44], the hypothesis of linearity is rejected in some traces (e.g. Tr5) and not in others (e.g. Tr3 and Tr4).

C. Katris, S. Daskalaki / Simulation Modelling Practice and Theory 75 (2017) 127–145

133

Table 3 Time dependent properties of real traces.

Trace

Stationaritya (ADF-test)

Tr1

Lag = 15 None: −5.6875 Intercept: −8.7382 Intercept+trend: −8.9462 (p-value < 0.01)

Tr2

Tr3

Tr4

Tr5

a

Lag = 15 None: −2.9672 Intercept: −8.4511 Intercept+trend: −8.8663 (p-value < 0.01) Lag = 15 None: −4.5573 Intercept: −9.8228 Intercept+trend: −9.9472 (p-value < 0.01) lag= 15 None: −3.0013 Intercept: −7.8074 Intercept + trend: −8.3516 (p-value< 0.01) Lag = 15 None: −4.0049 Intercept: −7.7504 Intercept+trend: −8.4254 (p-value< 0.01)

Autocorrelation (Ljung-Box)

LRD Hurst Exponent (H)

Non-linearity (White-NN test)

3496.5 (<0.01)

0.89767

7.2233 (0.07487)

2780 (<0.01)

0.82712

5.7184 (0.05731)

3026.6 (<0.01)

0.75789

3.0149 (0.2215)

3382.2 (<0.01)

0.86057

1.538 (0.4635)

3240.3 (<0.01)

0.82746

9.4833 (<0.01)

Stationarity is the alternative hypothesis.

Given that the data used for calculating the descriptive statistics are not i.i.d. and instead they are correlated, an effort has been made to estimate the bias introduced by the presence of autocorrelation using the block bootstrap approach [45]. This is a method for correcting estimations of moments when data are autocorrelated and stationary, like in our case. According to the block bootstrap the time series data yi are split in overlapping blocks of length b. Each block is represented as blockt = {yt , yt+1 , . . . , yt+b−1 }, where t = 1,…, N-b + 1 and N is the total number of observations of the time series, while resampling during bootstrap is realized using whole blocks. The crucial question for the block bootstrap technique is the decision about the length b of the blocks. In [46] some rules about the optimal block size are provided, however in this paper we borrow the concept of embedding dimension from dynamical systems. We consider a time series as dynamical system, where the value of each observation yi is a function of the m previous ones yt−1 , . . . , yt−m . Following the Embedding Theory of Takens and Mañé a sufficiently large embedding dimension m can capture the original time series and a widely used heuristic for calculating m is the False Nearest Neighbor (FNN) [47]. Detailed description for the time series as dynamical systems can be found in [48] and for their application on internet traffic forecasting in [30,49]. Specifically for the block bootstrap estimation technique the embedding dimension will be used as the size of the blocks. Using this method with 10,0 0 0 samples, we can estimate the bias in the time series and then subtract it from the initial biased estimations. Table 4 displays estimations for the mean, standard deviation, skewness and kurtosis both the initial (biased) and the block bootstrap ones, for all five traces examined previously. The last column gives the 90% confidence intervals calculated by the bootstrap method. We observe that for all cases examined the initial biased estimations lie inside the 90% bootstrap confidence intervals, suggesting that there are no significant differences between biased and corrected estimations. For this reason and due to the high computational effort required by the block bootstrap this step is skipped for the next steps.

5.2. Synthetic traffic using FARIMA models As a first step in creating and comparing traffic generators for the five traces just examined, we are fitting the FARIMA and FARIMA-GARCH models to the traffic traces. This task involves decisions on (i) the order of the corresponding ARIMA model, using BIC (Bayes Information criterion), (ii) the fractional parameter d, using the Geweke and Porter-Hudak (GPH) estimator, and (iii) the rest of the parameters of each model, using Maximum Likelihood. For the FARIMA-GARCH models, a GARCH (1,1) component is additionally appended to account for conditional heteroskedasticity. The fitting of the models has been realized with the help of R software [50,51] and for each video trace Table 5 displays the values adopted for the model parameters along with the BIC value for each one of the fitted model. For Tr4, the parameter d was set manually equal to 0.25 since the GPH estimator failed to indicate a value in the desired interval (0, 0.5).

134

C. Katris, S. Daskalaki / Simulation Modelling Practice and Theory 75 (2017) 127–145

Table 4 Block bootstrap estimations of moments and Hurst exponent. Silence of The Lambs (Embedding Dimension: 13) Moment Mean Sd Skewness Kurtosis

Biased 32.108 31.343 3.2556 19.438

Block bootstrap 31.396 29.071 2.942 17.130

90% Bootstrap CI 28.276- 34.232 23.849–33.92 2.103–3.591 9.250–23.815

Mr.Bean (Embedding Dimension: 5) Moment Biased Mean 31.608 Sd 16.9766 Skewness 1.922 Kurtosis 10.453

Block bootstrap 31.529 16.813 1.842 9.802

90% Bootstrap CI 30.026–33.055 15.097–18.51 1.317–2.350 6.058–13.459

Star Trek (Embedding Dimension: 8) Moment Biased Mean 17.640 Sd 12.799 Skewness 2.518 Kurtosis 14.158

Block bootstrap 17.736 12.738 2.436 13.520

90% Bootstrap CI 16.662–18.801 11.269–14.237 1.958–2.893 9.478–17.794

Die Hard III (Embedding Dimension: 7) Moment Biased Mean 42.528 Sd 22.150 Skewness 1.363 Kurtosis 5.325

Block bootstrap 42.759 22.055 1.354 5.293

90% Bootstrap CI 40.549–45.106 20.345–23.699 1.162–1.547 4.407–6.256

The Firm (Embedding Dimension: 7) Moment Biased Mean 14.309 Sd 10.915 Skewness 2.143 Kurtosis 10.112

Block bootstrap 14.299 10.829 2.107 9.838

90% Bootstrap CI 13.195–15.406 9.720–11.917 1.734–2.464 7.102–12.682

Table 5 Optimal parameter values for the models fitted to traces. Trace

Order (p, q)

Fractional parameter (d)

Tr1 FARIMA-N FARIMA– t FARIMA-GARCH-N FARIMA-GARCH – t

(3,1)

0.47113

Tr2 FARIMA-N FARIMA– t FARIMA-GARCH-N FARIMA-GARCH – t

(1,4)

Tr3 FARIMA-N FARIMA– t FARIMA-GARCH-N FARIMA-GARCH – t

(1,3)

Tr4 FARIMA-N FARIMA– t FARIMA-GARCH-N FARIMA-GARCH – t

(1,2)

Tr5 FARIMA-N FARIMA– t FARIMA-GARCH-N FARIMA-GARCH – t

(1,2)

a

BIC 7.629 6.837 6.970 6.503

0.32896 7.299 7.036 7.042 6.916 0.31316 6.487 6.165 6.216 6.039 0.2500a 7.159 6.725 6.945 6.631 0.49882

Parameter d was set manually equal to 0.25.

5.933 5.608 5.594 5.430

C. Katris, S. Daskalaki / Simulation Modelling Practice and Theory 75 (2017) 127–145

135

Table 6 Statistical characteristics of generated traffic traces. Trace

# of data

Mean

Sd

Skewness

Kurtosis

H

Ljung-Box test

White NN test

T1 FARIMA-N FARIMA– t FARIMA-GARCH–N FARIMA-GARCH – t

40 0 0 16,845 25,055 2597 24,559

32.108 22.820 33.801 29.440 48.710

31.343 17.605 33.565 31.813 60.164

3.256 1.158 3.381 2.561 3.425

19.438 4.242 26.634 10.001 17.260

0.898 0.785 0.816 0.760 0.975

3496.5 (<0.01) 11,717 (<0.01) 22,494 (<0.01) 1759.3 (<0.01) 24,013 (<0.01)

7.223 (0.075) 504.29 (<0.01) 315.42 (<0.01) 45.955 (<0.01) 96.053 (<0.01)

T2 FARIMA-N FARIMA– t FARIMA-GARCH–N FARIMA-GARCH – t

40 0 0 91,527 86,672 82,336 95,839

31.608 27.518 28.983 23.554 35.356

16.977 14.572 17.544 14.576 18.211

1.922 0.451 1.166 1.067 1.787

10.453 2.863 7.261 2.379 13.331

0.827 0.803 0.823 0.899 0.864

2780 (<0.01) 56,970 (<0.01) 63,984 (<0.01) 53,120 (<0.01) 75,307 (<0.01)

5.7184 (0.057) 400.06 (<0.01) 330.73 (<0.01) 481.4 (<0.01) 188.88 (<0.01)

T3 FARIMA-N FARIMA– t FARIMA-GARCH–N FARIMA-GARCH – t

40 0 0 83,481 75,075 81,216 87,876

17.640 17.002 19.128 15.344 19.689

12.799 10.172 14.356 9.978 13.415

2.518 0.619 2.861 1.278 3.831

14.158 3.063 31.820 3.840 58.231

0.758 0.770 0.780 0.811 0.848

3026.6 (<0.01) 55,104 (<0.01) 59,282 (<0.01) 54,165 (<0.01) 73,239 (<0.01)

3.0149 (0.221) 741.26 (<0.01) 373.42 (<0.01) 451.34 (<0.01) 269.4 (<0.01)

T4 FARIMA-N FARIMA – t FARIMA-GARCH-N FARIMA-GARCH – t

40 0 0 95,914 84,745 96,153 94,273

42.528 39.960 44.355 39.865 42.527

22.150 18.840 29.392 18.049 25.548

1.363 0.303 2.202 0.541 6.576

5.325 -0.221 15.225 0.915 139.406

0.861 0.762 0.785 0.767 0.822

3382.2 (<0.01) 76,093 (<0.01) 76,128 (<0.01) 74,702 (<0.01) 85,437 (<0.01)

1.538 (0.464) 178,79 (<0.01) 174.3 (<0.01) 217.92 (<0.01) 100.87 (<0.01)

T5 FARIMA-N FARIMA – t FARIMA-GARCH-N FARIMA-GARCH – t

40 0 0 11,286 12,319 18,447 88,520

14.309 8.425 13.010 9.232 34.873

10.915 6.476 12.499 8.199 20.564

2.143 1.184 2.253 2.146 1.031

10.112 1.328 8.170 7.724 4.287

0.827 0.805 0.854 0.849 0.973

3240.3 (<0.01) 6982.1 (<0.01) 9665.2 (<0.01) 12,525 (<0.01) 85,103 (<0.01)

9.4833 (<0.01) 303.03 (<0.01) 290.81 (<0.01) 215.51 (<0.01) 48.513 (<0.01)

Given the fitted models, we next simulate 10 0,0 0 0 values from each model, out of which we keep only those that are larger than the minimum value of each trace. The procedure is repeated four times for each trace, in order to generate from FARIMA and FARIMA-GARCH with Normal and Student’s t errors. In Table 6 we present the final size of the simulated traces, their descriptive statistics (mean, standard deviation, skewness and kurtosis), values of Hurst exponents (H) using the Aggregated Variance method, the results from the Ljung-Box tests for autocorrelation and White’s NN test for non-linearity. For readability purposes we also report the corresponding values for the real traces shown in Table 2. Using the BIC values in Table 5, one could conclude that FARIMA-GARCH with Student’s t distribution for the residuals gives the best model for all five datasets. However, in Table 6, this conclusion is overturned with the calculated descriptive statistics. It is further observed that the traces generated from the classical FARIMA models (FARIMA-N) underestimate the mean, variance, skewness and kurtosis, but capture reasonably well the estimations for H. On the other hand, the traces generated by FARIMA with Student’s t innovations (FARIMA-t) give descriptive statistics closer to those of real traces, even though sometimes excessively overestimated, especially with kurtosis. Traces generated from classical FARIMA-GARCH models (FARIMA-GARCH-N) in some cases display a slight improvement over the FARIMA-N in capturing descriptive statistics; however they still underestimated most parameters. Finally, FARIMA-GARCH with Student’s t innovations (FARIMA-GARCH-t) generated traces that significantly overestimated the descriptive statistics, for example the mean and sd for Silence of the Lambs (Tr1) or the kurtosis in Star Trek (Tr3) and Die Hard III (Tr4). We should also notice that all models generated traffic for which the assumption of linearity is rejected, while this was not the case for at least two of the traces. On the other hand, the Ljung-Box tests for autocorrelation applied on the synthetic traces agree with the corresponding results on the real traces. Judging from the statistics calculated for the synthetic traces, FARIMA with Student’s t innovations and FARIMA-GARCH with normal innovations appear to be closer to those calculated from the real traces. Therefore we expect that they would be more successful in generating synthetic traffic that mimic closely the behavior of each corresponding video. The next step was to verify this assertion by letting real and synthetic traces compete in the performance evaluation scenario presented in Section 4, compare them and then draw conclusions about the accuracy and appropriateness of the models. Each trace is evaluated based on its performance in queuing situations where either the buffer size is held constant at 100 (kbytes) and the service rate varies as percentage of a predefined reference value equal to the estimated mean of the trace plus one standard deviation or where the service rate is held constant at a level equal to the 80% of the reference value and the buffer size varied from 50 to 450 (kbytes). For both evaluations the performance measure is the loss probability measured from the simulation run with the corresponding queuing system. Fig. 1 displays the results of those simulations with the queuing systems, when each one of the five real traces was used as input process as well as each of the corresponding four synthetic traces generated from the different time series models. According to Fig. 1, the traces generated by the FARIMA-N models always underestimated the loss probability and captured the queuing behavior of the actual trace in very few situations, e.g. the Star Trek trace and only for high service rates.

136

C. Katris, S. Daskalaki / Simulation Modelling Practice and Theory 75 (2017) 127–145

Fig. 1. Performance evaluation of FARIMA and FARIMA-GARCH traffic generators.

C. Katris, S. Daskalaki / Simulation Modelling Practice and Theory 75 (2017) 127–145

137

Fig. 2. Empirical distributions for the real traces compared to Lognormal models.

The behavior of the traces generated by the FARIMA-GARCH-N models in most cases is close to that of the FARIMA-N models, and only in the case of Silence of the Lambs was very close to the actual trace. On the contrary, FARIMA-GARCH-t overestimated significantly the loss probability (except of the trace for Die Hard III) and so it appears to be unsuitable as a model, despite of the smallest BIC values it gave in fitting. The model with the best overall generators appears to be the FARIMAt, since it gives better estimates compared to the others. Still its performance is far from being considered as completely satisfying. Specifically, when buffer size varies loss probabilities almost always are overestimated with FARIMA-t, and when service rate varies then FARIMA-t follows the general trend of loss probabilities but only in few occasions it fully captures the queuing behavior of the real traffic. For this reason in the next section we attempt a new approach that involves first the generation of a data series by a FARIMA model and afterwards its transformation to lognormal variates through its inverse CDF function. 5.3. Synthetic traffic using the 2-step procedure The next stage of our analysis was to compare the “best” so far model, shown previously to be the FARIMA-t, with the 2step procedure for generating traffic presented in Section 3. Before doing so, however, the empirical distribution is calculated for all traces and compared to Lognormal as depicted in Fig. 2. For the construction of the empirical distributions there are a few different kernel functions that can be used, such as uniform, triangular, Epanenchnikov and Gaussian. For this work, the Gaussian kernels were used and for the bandwidth h the univariate plug-in selection method of [40], respectively. The fitting of the Lognormal models for each trace was realized using the PDF of LN (μ ˆ log yi , σˆ log yi ). It is quite clear from Fig. 2 that the Lognormal models for all five traces fit very satisfactorily on the empirical distributions plotted, so these were our choices as target distributions for the 2nd step of the procedure that follows. The procedure starts by generating a trace using a classical FARIMA model and afterwards all values are transformed into Lognormal variates, which form a new synthetic trace with different characteristics compared to the original one. We now want to compare the statistical characteristics and queuing performance of the new traces (named FARIMA-toLN) with the most successful model of the FARIMA-based models (FARIMA-t) and two benchmarks (ARIMA and Lognormal). We consider as benchmark models the ARIMA, which fails to capture the LRD property but it is considered as a very good choice for capturing SRD, and the lognormal distribution which generates independent data. For each model all parameters

138

C. Katris, S. Daskalaki / Simulation Modelling Practice and Theory 75 (2017) 127–145 Table 7 Statistical characteristics of synthetic traces. Trace

BIC

# of data

Mean

Sd

Skewness

Kurtosis

H

Ljung-box test

White NN test

T1 Lognormal ARIMA FARIMA- t FARIMA-to-LN

8.649 7.668 6.837 6.612

40 0 0 10 0,0 0 0 82,060 25,055 99,997

32.108 31.563 40.669 33.801 31.750

31.343 29.536 24.551 33.565 29.310

3.256 3.707 0.619 3.381 3.010

19.438 29.114 0.033 26.634 15.676

0.898 0.488 0.626 0.816 0.913

3496.5 (<0.01) 1.071 (0.301) 66,620 (<0.01) 22,494 (<0.01) 89,062 (<0.01)

7.223 (0.075) 4.367 (0.113) 680.44 (<0.01) 315.42 (<0.01) 688.34 (<0.01)

T2 Lognormal ARIMA FARIMA- t FARIMA-to-LN

8.191 7.304 7.036 5.440

40 0 0 10 0,0 0 0 96,031 86,672 10 0,0 0 0

31.608 31.746 32.590 28.983 31.813

16.977 17.929 15.472 17.544 17.879

1.922 1.925 0.295 1.166 1.895

10.453 7.276 -0.246 7.261 7.129

0.827 0.489 0.699 0.823 0.816

2780 (<0.01) 0.315 (0.574) 62,327 (<0.01) 63,984 (<0.01) 65,906 (<0.01)

5.7184 (0.057) 2.934 (0.231) 190.79 (<0.01) 330.73 (<0.01) 718.37 (<0.01)

T3 Lognormal ARIMA FARIMA- t FARIMA-to-LN

7.381 6.494 6.165 5.015

40 0 0 10 0,0 0 0 90,480 75,075 10 0,0 0 0

17.640 17.835 19.726 19.128 17.876

12.799 13.930 10.827 14.356 13.889

2.518 2.898 0.462 2.861 2.905

14.158 17.238 -0.147 31.820 18.605

0.758 0.489 0.616 0.780 0.790

3026.6 (<0.01) 0.736 (0.391) 62,065 (<0.01) 59,282 (<0.01) 70,896 (<0.01)

3.0149 (0.221) 3.789 (0.15) 447.12 (<0.01) 373.42 (<0.01) 925.39 (<0.01)

T4 Lognormal ARIMA FARIMA-t FARIMA-to-LN

8.745 7.169 6.725 5.281

40 0 0 10 0,0 0 0 96,449 84,745 10 0,0 0 0

42.528 42.626 43.670 44.355 42.705

22.150 23.415 20.314 29.392 23.279

1.363 1.862 0.269 2.202 1.8064

5.325 6.789 -0.253 15.225 6.330

0.861 0.489 0.666 0.785 0.769

3382.2 (<0.01) 0.290 (0.59) 79,090 (<0.01) 76,128 (<0.01) 79,992 (<0.01)

1.538 (0.464) 2.8783 (0.237) 160.21 (<0.01) 174.3 (<0.01) 432.14 (<0.01)

T5 Lognormal ARIMA FARIMA-t FARIMA-to-LN

6.954 5.949 5.608 4.328

40 0 0 10 0,0 0 0 89,059 12,319 10 0,0 0 0

14.309 14.291 16.310 13.010 14.414

10.915 11.375 9.072 12.499 11.596

2.143 2.972 0.491 2.253 2.541

10.112 18.184 -0.124 8.170 10.661

0.827 0.489 0.631 0.854 0.940

3240.3 (<0.01) 0.768 (0.381) 66,360 (<0.01) 9665.2 (<0.01) 88,568 (<0.01)

9.4833 3.8488 503.25 290.81 705.23

(<0.01) (0.146) (<0.01) (<0.01) (<0.01)

are determined with a typical fitting procedure to the actual traces and then from each model we simulate 10 0,0 0 0 values to form synthetic traffic. Table 7 displays the statistical characteristics of the synthetic traces generated by the different models. Specifically it shows the BIC values from the model fitting procedure, the final size of each trace after cleaning the simulated data, the descriptive statistics, H exponents, and outcomes of Ljung-Box and White NN tests for autocorrelation and nonlinearity. The BIC value for the FARIMA-to-LN models is calculated approximately as a function of sample size instead of the log-likelihood formula used in all other cases. We observe that the FARIMA-to-LN models always give the smallest BIC value, however since this is only an approximation we cannot extract any useful conclusion from it. Comparison of the statistical characteristics indicates that the traces generated from the FARIMA-to-LN models in most cases approximate fairly well the corresponding values of the actual traces, but so do also the traces generated from the FARIMA-t. An additional comparison is presented Fig. 3, where we compare again the four models for their ability to capture the queuing behavior of real traces across a range of buffer sizes and service rates. We observe that the synthetic traces generated from the procedure FARIMA-to-LN emulate closer the queuing behavior of the real traces by estimating quite accurately the loss probability in most cases. An exception occurs at The Firm trace, where FARIMA-to-LN overestimates the loss only across buffer sizes and then FARIMA-t is more accurate. This is not the case, however, for lower service rates where again the FARIMA-to-LN trace is closer to the real trace. In Fig. 3 we observe that overall the FARIMA-to-LN approximates the queuing behavior of the real traffic consistently and is more accurate than the other models in almost all cases. With Lognormal distribution, we take into account only the PDF of the data and we omit the dependence structure, i.e. SRD and LRD. In this case the model underestimates loss probability, and the accuracy is very low. FARIMA-t and ARIMA models only in rare occasions seem to approach the real traffic, e.g. ARIMA is close with Mr. Bean and Star Trek traces only for large buffer sizes and FARIMA-t is close to the Firm trace across service rates. 5.4. Experimental evaluation of dependence structure In this section we experiment further with the 2-step procedure applied previously. More specifically our interest is focused on whether the projection of a data series to a Lognormal distribution through its inverse CDF maintains, and at what level, the dependence structure that existed in the original series. For this reason, we consider simulated traces either completely random from a Normal(0, 1) distribution, or with some short range dependence structure from an ARIMA(1, 1) and an ARIMA(1, 1)-GARCH(1, 1) model, or with some long range dependence structure from a FARIMA(1, d, 1) and a FARIMA(1, d, 1)-GARCH(1, 1)model. For each one of these models we simulated eight traces of different sizes and then applied the transformation to the Lognormal with parameters 0 and 1. The sizes of the traces varied from 10 0 0 to 50 0,0 0 0 observations. The means and standard deviations are calculated for all transformed traces and are depicted in Fig. 4. It is obvious that samples are required to be of size approximately 10 0,0 0 0 in order for the first two moments to converge to certain

C. Katris, S. Daskalaki / Simulation Modelling Practice and Theory 75 (2017) 127–145

139

Silence of the Lambs (Tr1) 0,6

0,3 0,25

0,4

Loss

Loss

0,2 0,15

0,2

0,1 0,05

0

0 50

100

150

200

250

300

350

400

0,4

450

0,5

0,6

0,7

0,8

0,9

1

0,8

0,9

1

Service Rates

Buffer Size

Mr. Bean (Tr2) 0,8

0,30 0,25

0,6

Loss

Loss

0,20 0,15 0,10

0,4

0,2

0,05 0

0,00 50

100

150

200

250

300

350

400

0,4

450

0,5

0,6

0,7 Service Rates

Buffer Size

Star Trek (Tr3) 0,6

0,30 0,25

0,4

Loss

Loss

0,20 0,15

0,2

0,10 0,05

0

0,00 50

100

150

200

250

300

350

400

0,4

450

0,5

0,6

0,7

0,8

0,9

1

Service Rates

Buffer Size

Die Hard III (Tr4) 0,8

0,30 0,25

0,6

Loss

Loss

0,20 0,15 0,10

0,4

0,2

0,05 0

0,00 50

100

150

200

250

300

350

400

0,4

450

0,5

0,6

0,7

0,8

0,9

1

0,8

0,9

1

Service Rates

Buffer Size

The Firm (Tr5) 0,30

0,6

0,25 0,4

Loss

Loss

0,20 0,15

0,2

0,10 0,05

0

0,00 50

100

150

200

250 Buffer Size

300

350

400

450

0,4

0,5

0,6

0,7 Service Rates

Fig. 3. Performance evaluation of FARIMA-to-LN, benchmark models and FARIMA-t.

140

C. Katris, S. Daskalaki / Simulation Modelling Practice and Theory 75 (2017) 127–145

1,68

3

1,67

2,5 Standard Deviation

1,66

Mean

1,65 1,64 1,63 1,62 1,61 1,6

2 Normal

1,5

ARIMA

1

ARIMA/GARCH FARIMA

0,5

1,59

FARIMA/GARCH

1,58

0

Observations

12

600

10

500

8

400 Kurtosis

Skewness

Observations

6

Normal

300

ARIMA

4

200

ARIMA/GARCH

2

100

FARIMA FARIMA/GARCH

0

0

observations

observations Fig. 4. Moments of the transformed traces.

values. On the contrary, skewness and kurtosis, require samples of larger size in order to converge to some value, especially for the FARIMA-GARCH model. In addition, the Hurst exponents are calculated for the original samples and for the transformed ones and their values, along with the % change between them, are reported in Table 8 and depicted in Fig. 5. One may observe that the random generator for Normal variates produces a sample of independent numbers, so the H exponent is approximately 0.5 both for the original and the transformed traces. Similarly, the ARIMA (1,1) and ARIMA-GARCH (1,1)/(1,1) simulators produce traces that do not carry long memory and H converges to some value around 0.51, while the transformed traces follow very close the original values. Lastly, the FARIMA (1,0.25,1) and FARIMA-GARCH (1,0.25,1) /(1,1) simulators produce traces that display LRD with the H exponent to converge to some value close to 0.75. When the transformation to lognormal distribution is applied the new traces lose some of it and converge to a smaller value, around 0.71 and 0.72, respectively, however still indicating LRD. Looking at Fig. 5(a) and (b) it is obvious that the three first simulators converge to their corresponding H values with samples 10,0 0 0 observations and above, however for FARIMA and FARIMA/GARCH it requires at least 50,0 0 0. The detailed values for the aforementioned comparison are given in Table 8, where the percent differences between the H estimates for the original and transformed data series are also given, where Fig. 5(c) depicts the corresponding absolute percent differences. The experimentation just performed indicates that the transformed traces do maintain a major part of the dependence structure they carried prior to their projection. Apparently in order to generate traces that mimic closely real traffic it is necessary to start with a trace that carries the desired dependencies, measured by H, and the transformation through the inverse CDF of a suitable Lognormal distribution can do the rest. Another concern that this experimentation raised was the characterization of the difference observed between the H estimates for the original and the transformed data series. For this reason, and since it is known that there is not a single “best” technique for estimating the Hurst exponent, we extended the investigation to different techniques for estimating it. More specifically, we estimated H for the actual and the synthetic traces generated from the 2-step procedure using 6 different techniques: R/S, Aggregate Variance, Periodogram, Higutchi and Discrete Wavelet transform. Table 9 displays the estimation for the actual traces, for the simulated traces with FARIMA-to-LN and their absolute difference for each one technique, while Table 10 displays the average absolute differences between actual and simulated traces recorded for each estimation technique. It is observed that the average absolute differences between the estimations of H for the original and the transformed traces vary between 0.0027 (Higuchi method) and 0.0919 (R/S) for the different techniques. Using the Higuchi method the absolute difference never exceeded 0.009, using the Wavelet transforms method never exceeded 0.11,

C. Katris, S. Daskalaki / Simulation Modelling Practice and Theory 75 (2017) 127–145

141

Table 8 Statistical characteristics of original and transformed traces. Size of dataset Normal (0,1)

H (original)

H (transformed)

% change

10 0 0 10,0 0 0 50,0 0 0 10 0,0 0 0 20 0,0 0 0 30 0,0 0 0 40 0,0 0 0 50 0,0 0 0

0.54475 0.49554 0.48607 0.49257 0.49555 0.49547 0.49532 0.49599

0.62510 0.48491 0.50024 0.49502 0.49589 0.49710 0.49626 0.49775

14.75% -2.15% 2.91% 0.50% 0.07% 0.33% 0.19% 0.35%

0.53353 0.48902 0.49100 0.50228 0.50807 0.51003 0.50856 0.51016

14.08% −0.73% −2.50% −1.37% −1.09% −0.37% −0.60% −0.38%

ARIMA (1,1) – ar1 = 0.3, ma1 = 0.2, 10 0 0 10,0 0 0 50,0 0 0 10 0,0 0 0 20 0,0 0 0 30 0,0 0 0 40 0,0 0 0 50 0,0 0 0

μ = 0, σ = 1

ARIMA/GARCH (1,1) / (1,1) -μ = 0, 10 0 0 10,0 0 0 50,0 0 0 10 0,0 0 0 20 0,0 0 0 30 0,0 0 0 40 0,0 0 0 50 0,0 0 0

ar1 = 0.3, ma1 = 0.2, ω = 0.0 0 0 01, 0.47748 0.49210 0.50371 0.50946 0.51388 0.51187 0.51169 0.51235

0.46769 0.49261 0.50357 0.5092 0.51369 0.51194 0.51162 0.51209

FARIMA (1,d,1) - ar1 = 0.3, ma1 = 0.2, d = 0.25, μ = 0, σ = 1 10 0 0 0.73209 10,0 0 0 0.71606 50,0 0 0 0.73352 10 0,0 0 0 0.74469 20 0,0 0 0 0.74043 30 0,0 0 0 0.74923 40 0,0 0 0 0.75011 50 0,0 0 0 0.75097 FARIMA/GARCH (1,d,1) / (1,1) - ar1 = 0.3, ma1 = 0.2, μ = 0, α 1 = 0.05, 10 0 0 0.72709 10,0 0 0 0.71483 50,0 0 0 0.73369 10 0,0 0 0 0.74511 20 0,0 0 0 0.74082 30 0,0 0 0 0.74943 40 0,0 0 0 0.75020 50 0,0 0 0 0.75114

α 1 = 0.05, β 1 = 0.05 0.53634 0.49384 0.49292 0.50448 0.50870 0.51046 0.50919 0.51112

12.33% 0.35% −2.14% −0.98% −1.01% −0.28% −0.49% −0.24%

0.71466 0.66313 0.69505 0.70963 0.70816 0.71526 0.71384 0.71606

−2.38% −7.39% −5.25% −4.71% −4.36% −4.53% −4.83% −4.65%

β 1 = 0.05, ω = 0.00001, d = 0.25 0.71313 0.65209 0.68631 0.69811 0.69879 0.70535 0.70449 0.70612

−1.92% −8.78% −6.46% −6.31% −5.67% −5.88% −6.09% −5.99%

Table 9 Hurst Exponent estimations with different methods (actual traces, simulated traces, absolute difference). Trace

Wavelet

R/S

Aggregated variance

Periodogram

Higuchi

Average

Silence of The Lambs

1.0473 1.0475 0.0 0 02

0.9173 0.8424 0.0749

0.8977 0.9125 0.0148

1.0602 0.9962 0.0640

0.9602 0.9692 0.0090

0.97654 0.95356 0.0326

0.8549 0.9362 0.0813

0.7148 0.7938 0.0790

0.8271 0.8164 0.0107

0.9088 0.9745 0.0657

0.9702 0.9690 0.0012

0.85516 0.89798 0.0476

0.9860 0.9100 0.0760

0.8156 0.7822 0.0334

0.7579 0.7905 0.0326

0.9783 0.8571 0.1212

0.9704 0.9695 0.0 0 09

0.90164 0.86186 0.0528

1.1142 1.0069 0.1073

1.0 0 06 0.8530 0.1476

0.8606 0.7692 0.0914

1.0193 0.9631 0.0562

0.9700 0.9692 0.0 0 08

0.99294 0.91228 0.0807

1.1127 1.0360 0.0767

0.7049 0.8297 0.1248

0.8275 0.9398 0.1123

1.0222 1.0055 0.0167

0.9671 0.9688 0.0017

0.92688 0.95596 0.0664

Abs. Diff. Mr. Bean Abs. Diff. Star Trek Abs. Diff. Die Hard III Abs. Diff. The Firm Abs. Diff.

142

C. Katris, S. Daskalaki / Simulation Modelling Practice and Theory 75 (2017) 127–145

(b)

0,8

0,8

0,7

0,7

0,6

0,6

H (transformed)

H (original)

(a)

0,5 0,4 0,3 0,2

0,5 0,4 0,3 0,2

0,1

0,1

0

0

Observaons

Observaons

(c) 16,00%

H (% change)

14,00% 12,00% 10,00%

Normal

8,00%

ARIMA

6,00%

ARIMA/GARCH

4,00% 2,00%

FARIMA

0,00%

FARIMA/GARCH

Observaons Fig. 5. Hurst exponent values for (a) initial and (b) transformed data; (c) absolute % change in H. Table 10 Absolute Differences per method. Method

Average Absolute Difference

Wavelet R/S Aggregated Variance Periodogram Higuchi Average of all Methods

0.0683 0.0919 0.0524 0.0648 0.0027 0.0560

while the worst case was 0.15 (with R/S). The average of all absolute differences was found to be 0.056, so we conclude that the LRD measured in the simulated traces was very close to that of the original. 5.5. Generating traffic for a new video The methodology that was described in Section 3 and used for the experimentation presented previously has been shown to capture satisfactorily the static and time-dependent statistical characteristics of the corresponding real traces, resulting eventually to similar queuing behaviors when fed to a FIFO queuing system. The challenge now is to use the methodology in more practical situations where the real trace is not available and instead there is only an estimate (or a target value) for the mean and standard deviation. Since it is not necessary for the procedure to track down the exact LRD structure in a video trace in order to be successful then adopting a “typical” FARIMA model as reference would probably be sufficient. As discussed earlier, in our case we adopted the FARIMA (1,0.3,1) as generator for the initial trace and then projected it to

C. Katris, S. Daskalaki / Simulation Modelling Practice and Theory 75 (2017) 127–145

Fig. 6. Queuing performance of FARIMA(r)-to-LN.

143

144

C. Katris, S. Daskalaki / Simulation Modelling Practice and Theory 75 (2017) 127–145 Table 11 Statistical characteristics of the traces generated by the FARIMA(r)-to-LN procedure. Trace

Mean

Sd

Skewness

Kurtosis

H

Ljung-Box test

White NN test

T1 FARIMA(r)-to-LN Abs. Difference

32.108 32.138 0.030

31.343 31.663 0.320

3.256 4.147 0.891

19.438 44.568 25.130

0.898 0.771 0.127

3496.5 (<0.01) 67,914 (<0.01)

7.223 (0.075) 1296.2 (<0.01)

T2 FARIMA(r)-to-LN Abs. Difference

31.608 31.610 0.002

16.977 17.023 0.046

1.922 1.822 0.100

10.453 9.678 0.775

0.827 0.786 0.041

2780 (<0.01) 72,592 (<0.01)

5.7184 (0.057) 594.64 (<0.01)

T3 FARIMA(r)-to-LN Abs. Difference

17.640 17.644 0.004

12.799 12.867 0.068

2.518 2.678 0.160

14.158 18.582 4.424

0.758 0.780 0.022

3026.6 (<0.01) 70,896 (<0.01)

3.0149 (0.221) 925.39 (<0.01)

T4 FARIMA(r)-to-LN Abs. Difference

42.528 42.531 0.003

22.150 22.207 0.057

1.363 1.756 0.393

5.325 9.166 3.841

0.861 0.786 0.074

3382.2 (<0.01) 72,734 (<0.01)

1.538 (0.464) 568.35 (<0.01)

T5 FARIMA(r)-to-LN Abs. Difference

14.309 14.312 0.003

10.915 10.979 0.064

2.143 2.870 0.727

10.112 21.192 11.080

0.827 0.778 0.049

3240.3 (<0.01) 70,361 (<0.01)

9.4833 (<0.01) 924.2 (<0.01)

a Lognormal distribution with parameters that are estimated using the MCMC routine from the bayescount R package [52]. Using the five traces T1-T5 we apply this idea and compare the performance of the newly generated traces (FARIMA(r)-to-LN) with the actual ones. Table 11 displays the descriptive statistics, Hurst exponents for the traces generated by the FARIMA(r)-to-LN procedure and for the real traces for comparison; in addition, the outcomes from the Ljung-Box and White NN tests. With the exception of kurtosis, in most cases the descriptive statistics calculated for the synthetic traces are very close to the corresponding values calculated for the real traces. The H exponent in most cases is different than the values we saw in Table 7, since now the FARIMA model is just a reference model, however it still indicates significant LRD for all traces. Lastly, the outcomes of the Ljung-Box and White NN tests for the simulated traces indicate autocorrelation and non-linear dependencies for all traces. With such proximity in most values of the statistical parameters observed in Table 11 it is expected that the synthetic traces will emulate closely the real traces in a queuing situation. This expectation is verified in Fig. 6 where the synthetic traces FARIMA(r)-to-LN are fed as input to the FIFO queuing system described in Section 4. Again the performance measure is the loss probability and is compared with the performance of the real traces. Their absolute difference in estimating loss probability never exceeded the value of 0.06, indicating that the two traces have very similar performance. 6. Summary and conclusions In this paper a 2-step procedure based on time series and the Inverse Sampling method is proposed in order to generate synthetic video traffic. The basis for the procedure is the classical FARIMA model and the simulated traces are transformed using the inverse CDF of a “target” marginal distribution, Lognormal in our case. The result is synthetic traffic with features such as autocorrelation, SRD, LRD and marginal distribution that is close to distributions observed in real video traffic. The procedure has been applied to five different publicly available video traces and the synthetic traces produced displayed behavior closer to that of the real traces, when compared to other generators in the FARIMA and FARIMA-GARCH family with normal or Student’s t innovations, or to other benchmarking models such as the and the Lognormal distribution. Finally, we considered the case of a new video trace for which we may only have estimation or a guess for its mean and standard deviation and not the actual trace. In such a case the proposed procedure is adjusted to start with a typical FARIMA model and then be projected to the desired Lognormal distribution. Apparently such an adjustment does not nullify the usefulness of the suggested procedure, on the contrary it is shown experimentally that the synergy between the FARIMA reference model and the Lognormal distribution can provide traces that emulate quite satisfactorily the behavior of the real traces. References ´ [1] E. Casilari, A. Reyes, A. Diaz-Estrella , F. Sandoval, Classification and comparison of modelling strategies for VBR video traffic, in: Proceedings of International Teletraffic Congress (ITC-16)’99, 1999. [2] A. Alheraish, Autoregressive video conference models, Int. J. Netw. Mgmt. 14 (5) (2004) 329–337. ´ [3] E. Casilari, A. Reyes, A. Diaz-Estrella , F. Sandoval, Characterisation and modelling of VBR video traffic, Electron. Lett. 34 (10) (1998) 968–969. [4] S. Tanwir, H. Perros, A survey of VBR video traffic models, IEEE Commun. Surv. Tutorials 15 (4) (2013) 1778–1802. [5] M. Krunz, S.K. Tripathy, On the characterization of VBR MPEG streams, Perform. Eval. Rev. 25 (1) (1997) 192–202. [6] H. Koumaras, C. Skianis, A. Kourtis, Analysis and modeling of H.264 unconstrained VBR video traffic, Int. J. Mob. Comput. Multimedia Commun. 1 (4) (2009) 14–31 IGI-Global. [7] D. Liu, E.I. Sára, W. Sun, Nested auto-regressive processes for MPEG-encoded video traffic modeling. IEEE Transactions on Circuits and Systems for Video Technology, 11(2), 169-183. [8] D.P. Heyman, T.V. Lakshman, Source models for VBR broadcast-video traffic, IEEE/ACM Trans. Netw. 4 (1) (1996) 40–48. [9] A. Lazaris, P. Koutsakis, Modeling multiplexed traffic from H.264/AVC videoconference streams, Comput. Commun. 33 (10) (2010) 1235–1242.

C. Katris, S. Daskalaki / Simulation Modelling Practice and Theory 75 (2017) 127–145

145

[10] M. Frey, S. Nguyen-Quang, A gamma-based framework for modeling variable-rate MPEG video sources: the GOP GBAR model, IEEE/ACM Trans. Netw. 8 (6) (20 0 0) 710–719. [11] Q.T. Zhang, A general AR-based technique for the generation of arbitrary gamma VBR video traffic in ATM networks, IEEE Trans. Circuits Syst. Video Technol. 9 (7) (1999) 1130–1137. [12] A. Alheraish, S.A. Alshebeili, T. Alamri, A GACS modeling approach for MPEG broadcast video, IEEE Trans. Broadcast. 50 (2) (2004) 132–141. [13] J.-L.C. Wu, Y.-W. Chen, C.-C. Shiu, Traffic modeling and bandwidth allocation for MPEG video sources in ATM networks, in: IEEE Global Telecommunications Conference, 3, 1995, pp. 2237–2241. [14] M.W. Garrett, W. Willinger, Analysis, modeling and generation of self-similar VBR video traffic, in: ACM SIGCOMM, 1994, pp. 269–280. [15] C. Huang, M. Devetsikiotis, I. Lambadaris, A.R. Kaye, Modeling and simulation of self-similar variable rate compressed video: a unified approach, in: ACM SIGCOMM, 1995, pp. 114–125. [16] H. Liu, N. Ansari, Y.Q. Shi, Modeling VBR video traffic by Markov-modulated self-similar processes, in: IEEE 3rd Workshop on Multimedia Signal Processing, 1999, pp. 363–368. [17] N. Ansari, H. Liu, Y.Q. Shi, H. Zhao, On modeling MPEG video traffics, IEEE Trans. Broadcast. 48 (4) (2002) 337–347. [18] U.K. Sarkar, S. Ramakrishnan, D. Sarkar, Modeling full-length video using Markov-modulated Gamma-based framework, IEEE/ACM Trans. Netw. 11 (4) (2003) 638–649. [19] Y. Sun, J.N. Daigle, A source model of video traffic based on full-length VBR MPEG4 video traces, in: GLOBECOM-IEEE Global Telecommunications Conference, 2, 2005, pp. 766–770. [20] M. Dai, Y. Zhang, D. Loquinov, A unified traffic model for MPEG-4 and H.264 video traces, IEEE Trancs. Multimedia 11 (5) (2009) 1010–1023. [21] M. Sheng, C. Ji, Modeling heterogeneous network traffic in wavelet domain, IEEE/ACM Trans. Netw. 9 (5) (2001) 634–649. [22] M. Cario, B. Nelson, Autoregressive to anything: time-series input processes for simulation, Oper. Res. Lett. 19 (1996) 51–58. [23] B. Biller, B. Nelson, Fitting time-series input processes for simulation, Oper. Res. 53 (2005) 549–559. [24] C. Granger, R. Joyeux, An introduction to long-memory time series models and fractional differencing, J. Time Ser. Anal. 1 (1980) 15–29. [25] J. Hosking, Fractional differencing, Biometrika 68 (1981) 165–176. [26] M. Corradi, R. Garroppo, S. Giordano, M. Pagano, Analysis of F-ARIMA processes in the modelling of broadband traffic, in: International Conference on Communications (ICC 2001), 3, IEEE, 2001, pp. 964–968. [27] R. Garroppo, S. Giordano, S. Lucetti, F. Russo, Comparison of LRD and SRD traffic models for the performance evaluation of finite buffer systems, in: International Conference on Communications, 2001 (ICC 2001), 9, IEEE, 2001, pp. 2681–2686. [28] G. Chiruvolu, R. Shankar, An approach towards resource management and transportation of VBR video traffic, in: International Conference on Communications (ICC ’97), IEEE, 1997, pp. 550–554. [29] N. Sadek, A. Khotanzad, T. Chen, ATM dynamic bandwidth allocation using F-ARIMA prediction model, in: International Conference on Computer Communications and Networks (ICCCN), IEEE, 2003, pp. 359–363. [30] C. Katris, S. Daskalaki, Combining time series forecasting methods for internet traffic, in: Stochastic Models, Statistics and their Applications, 2015, pp. 309–317. [31] J. Geweke, S. Porter-Hudak, The estimation and application of long memory time series models, J. Time Ser. Anal. 4 (1983) 221–238. [32] R. Baillie, C. Chung, M. Tieslau, Analysing inflation by the fractionally integrated ARFIMA-GARCH model, J. Appl. Econ. 11 (1996) 23–40. [33] S. Kim, Forecasting internet traffic by using seasonal GARCH models, J. Commun. Netw. 13 (2011) 621–624. [34] B. Zhou, D. He, Z. Sun, W. Ng, Network traffic modeling and prediction with ARIMA/GARCH, in: Proceedings Of HET-Nets Conference, 2005, pp. 1–10. [35] A. Syed, H. Saleem, H. Syed, MCMC simulation of GARCH model to forecast network traffic load, Int. J. Comput. Sci. Inf. 9 (2012) 277–284. [36] R. Engle, T. Bollerslev, Modelling the persistence of conditional variances, Econometric Rev. 5 (1986) 1–50. [37] I. Antoniou, V. Ivanov, V. Ivanov, P. Zrelov, On the log-normal distribution of network traffic, Physica D 167 (2002) 72–85. [38] B. Silverman, Density Estimation for Statistics and Data Analysis, Chapman and Hall, London, 1986. [39] T. Duong, ks: Kernel Smoothing, 2015. [40] M.P. Wand, C. Jones, Multivariate plug-in bandwidth selection, Comput. Stat. 9 (2) (1994) 97–116. [41] M.S. Taqqu, V. Teverovsky, W. Willinger, Estimators for long-range dependence: an empirical study, Fractals 3 (4) (1995) 785–798. [42] P. Abry, D. Veitch, Wavelet analysis of long-range-dependent traffic, IEEE trans. Inf. Theory 44 (1) (1998) 2–15. [43] F. Fitzek, M. Reisslein, MPEG-4 and H.263 video traces for network performance evaluation, IEEE Netw. 15 (2001) 40–54. [44] T. Lee, H. White, C. Granger, Testing for neglected nonlinearity in time series models, J. Econometrics 56 (1993) 269–290. [45] H. Kunsch, The jackknife and the bootstrap for general stationary observations, Ann. Statist. 17 (1989) 1217–1241. [46] P. Hall, J. Horowitz, B. Jing, On blocking rules for the bootstrap with dependent data, Biometrika 82 (1995) 561–574. [47] H. Abarbanel, M. Kennel, Local false nearest neighbors and dynamical dimensions from observed chaotic data, Phys. Rev. E 47 (1993) 3057–3068. [48] R. Frank, N. Davey, S. Hunt, Time series prediction and neural networks, J. Intell. Rob. Syst. 31 (2001) 91–103. [49] C. Katris, S. Daskalaki, Comparing forecasting approaches for Internet traffic, Expert Syst. Appl. 42 (2015) 8172–8183. [50] A. Ghalanos, rugarch: univariate GARCH models, 2015. [51] C. Fraley, F. Leisch, M. Maechler, V. Reisen, A. Lemonte, fracdiff: fractionally differenced ARIMA aka ARFIMA(p,d,q) models, 2012. [52] M. Denwood, Bayescount: power calculations and bayesian analysis of count distributions and fecrt data using mcmc, 2015