Efﬁcient Generators of the Generalized Fractional Gaussian Noise and Cauchy Processes

: In the last years of the past century, complex correlation structures were empirically observed, both in aggregated and individual trafﬁc traces, including long-range dependence, large-timescale self-similarity and multi-fractality. The use of stochastic processes consistent with these properties has opened new research ﬁelds in network performance analysis and in simulation studies, where the efﬁcient synthetic generation of samples is one of the main topics. Nowadays, networks have to support data services for trafﬁc sources that are poorly understood or still insufﬁciently observed, for which simple, reproducible, and good trafﬁc models are yet to be identiﬁed, and it is reasonable to expect that previous generators could be useful. For this reason, as a continuation of our previous work, in this paper, we describe efﬁcient and online generators of the correlation structures of the generalized fractional noise process (gfGn) and the generalized Cauchy (gC) process, proposed recently. Moreover, we explain how we can use the Whittle estimator in order to choose the parameters of each process that give rise to a better adjustment of the empirical traces.


Introduction
Fundamental network algorithms and key performance metrics in telecommunication networks and services, such as routing, delay, age of information, or buffer sizing, rely on accurate statistical traffic models capable of replicating the temporal and spatial correlation observable in many diverse packet streams [1,2]. Further, current networks have been architected to support data services for traffic sources that are poorly understood or still insufficiently observed, such as mMTC (massive machine-type communication) or URLLC (ultra reliable low-latency communications) in 5G/6G, for which simple, reproducible, and good traffic models are yet to be developed. Since, in the past, complex correlation structures were empirically observed both in aggregated and individual traffic traces, including long-range dependence (LRD) and large-timescale self-similarity (SS) [3][4][5][6][7][8][9][10][11][12][13][14], it is reasonable to expect that novel or modified flexible stochastic processes have to be analyzed for teletraffic and simulation analysis of such advanced communication services. A similar research effort was made when queuing models with fractal correlated input and the analysis of that on network performance [15][16][17][18][19][20][21][22][23][24][25][26] were thoroughly studied. These works demonstrated that second-order statistics, and their extreme manifestations as LRD and SS, lead to very slow decay of the queue backlogs, thus having a huge impact on loss and delay.
Apart from the specific application in communications networks, long-range dependent and self-similar processes appear in many other fields in science and engineering. In [27], the authors show examples in physics, chemistry and biology. They are also of interest in the analysis of several types of time series, such as meteorological and hydrological data, stock markets data, rare events in finance and insurance, or biometric signals [28] and recently in vehicular traffic data [29].
There are, in the literature, several classes of stochastic processes with LRD and/or SS properties which have been successfully used to model the complex time series of real traffic. These include fractional Gaussian noise [15] (fGn) and its recent generalization [30]; alpha-stable processes [31]; the Cauchy process and its generalization [32]; fractional ARIMA [33] (fARIMA); the sum of on-off sources with sojourn times following a heavytailed distribution [7]; the state process of a M/G/∞ queuing system [34]; and wavelet functions [35].
However, in addition to the analytical tractability and flexibility to give rise to different correlation functions, an important feature with these advanced models is computational efficiency in the numerical simulation of their sample paths, especially because when LRD or SS are present, very long synthetic sample paths are necessary for accurate simulations. In this regard, one major drawback of many direct generators is that these are efficient only because of the application of an off-line algorithm for the creation of the traces.
To overcome this drawback, in our previous work [36,37], we leveraged the singular properties of the occupancy process in the M/G/∞ queuing system to generate, at low computational cost, arbitrarily long sequences with arbitrary covariance function. The key idea is just to select the service time distribution of customers in the M/G/∞ according to a suitable discrete probability mass function, and take advantage of the memoryless property in the arrivals to implement a very fast discrete-event simulation of this elementary model. In this paper, we show how to select the service time distribution of the M/G/∞ system in order to obtain accurate and efficient generators of the generalized fGn process [38] and the generalized Cauchy process [39], two classes of processes that provide a powerful way to describe the multi-fractal phenomena of traffic [40] because they are indexed by two independent parameters.
As an application, we use these efficient generators to reproduce the covariance structure of a real traffic sequence. In order to estimate the parameters of the models, which are the input that the simulator needs, we apply the Whittle estimator [41][42][43], a sufficient and robust statistic that operates in the spectral domain.
The remainder of the paper is organized as follows. In Section 2, we review the main concepts related to short-range dependence (SRD), LRD and self-similarity. The generalized fGn process is presented in Section 3.1 and the generalized Cauchy process is Section 3.2. The main properties of the M/G/∞ process used in this work are described in Section 4. In Sections 5.1 and 5.2, we describe the improved M/G/∞-based generator of the correlation structure of the gfGn and generalized Cauchy processes, respectively and in Section 5.3 we show some results related to the accuracy of the generators. In Section 6, we explain how we can use the Whittle estimator for fitting empirical traffic. Finally, in Section 7, we summarize the conclusions.

SRD, LRD and Self-Similarity
We briefly review in this section the definitions and main properties of self-similar stochastic processes and long-range statistical correlation [44].
A stationary stochastic process X = (X n ) n≥0 is said to have short-range dependence (SRD) whenever the sum of the covariance coefficients is convergent, i.e., . A class of processes having SRD is that wherein the covariance series decays exponentially: An equivalent condition to exponential decay is that the spectral density f X (ω) := ∑ ∞ k=0 r k e −ωk is bounded at the origin. If the sum of the covariance series is not conver-gent, the process X is termed long-range dependent (LRD). An example of LRD processes is the broad class of processes having sub-exponential (e.g., hyperbolic) decay: so r k = o(k −α ). In this case, the spectral density of X has a singularity at the origin. Self-similarity is defined based on the properties of the covariance function of the aggregated process X (m) = (X (m) j ) j≥1 built from X when X has finite variance. This is the average of X over non-overlapping blocks of length m Then, X is called exactly second-order self-similar with similarity parameter H [45] if m 1−H X (m) has the same covariance function as X for all m. This means that the aggregated process shows the same statistical correlation over time as the original process X. The covariance function of a self-similar process is Thus, that is, r H k decays hyperbolically (as in (2)) with α = 2H − 2. Clearly, X possesses LRD when H ∈ (0.5, 1). The covariance function (4) implies a spectral density given by where c is a normalization constant such that π −π f H (ω)dω = Var[X]. When (4) holds asymptotically, and, in addition, the convergence condition lim m→∞ = r (m) k = r H k for all k ≥ 1 is satisfied, then the process is called asymptotically second-order self-similar. Clearly, for asymptotic self-similarity, the spectral density converges to that of a self-similar process lim m→∞ f X (m) (ω) = f H (ω) for any ω ∈ [−π, π]. It is not hard to see that a stationary process with a covariance function with hyperbolic decay, such as (2), is asymptotically second-order self-similar.
Stationary Gaussian processes are not always exactly or asymptotically self-similar, but it is known that they satisfy a property of local self-similarity as defined in [46]. The covariance function function r k is for k → 0, where 0 < α < 2 is the order of the process. The so-called fractal dimension is δ = 2 − α 2 . Therefore, the behavior of the covariance function near the origin determines the local irregularity of the process, and larger values of δ imply more irregularity (or burstiness) at small timescales.

Generalized Fractional Gaussian Noise
Let B H (t) denote a fractional Brownian motion process [47]. Then, its sequence of discrete increments X n := B H (n) − B H (n − 1), n = 1, 2, . . . is an exactly self-similar stationary Gaussian process, known as fractional Gaussian noise (fGn). Since self-similarity is exact, its covariance function and spectral density are given by (4) and (6), respectively. Despite its fixed covariance structure, the fGn process is known to be a good approximation of the aggregation of more complex LRD Gaussian and non-Gaussian processes [48,49].
In [38], the generalized fractional Gaussian noise (gfGn) is introduced, with parameters 0 < H < 1 and 0 < a ≤ 1. This is a Gaussian process with covariance function that can be approximated by r H,a k ≈ H(2H − 1)|k a | 2H−2 , so the LRD condition is clearly satisfied if 0.5 < H < 1. The spectral density function of gfGn is [50] f H,a (ω) = sin(Haπ)Γ(2Ha which can be approximated as The latter implies that f H,a (0) → ∞, which is actually the LRD condition reflected in the spectral domain. Figure 1 shows the analytical covariance function for different values of the parameters. It is increasing in H (for fixed a), and decreasing in a (for fixed H).

Generalized Cauchy Process
The generalized Cauchy process X C = (X n ) n≥1 with parameters 0 < α ≤ 2 and β > 0 is a stationary Gaussian process with an autocorrelation function given by [39,51] ρ β,α The covariance function ρ β,α k is positive-definite for the above ranges of α and β, and it is completely monotone for 0 < α ≤ 1 and β > 0. When α = β = 2, one obtains the usual Cauchy process. Since lim k→∞ k β ρ β,α Its spectral density function is [52] Here, and Letting ω ≈ 0, the behavior of the spectral density function (12) is well approximated by so f β,α (0) → ∞. This is, again, the LRD condition reflected in the spectral domain.

The M/G/∞ Process
The M/G/∞ process [53] is the stationary limit of the state of a M/G/∞ queuing system, a queuing model with Poisson arrivals, general independent and identically distributed service times distributed as the random variable S, and an unlimited number of servers. Instead of the usual view as a continuous-time process, it is more convenient in our case to discretize time and look at the number of servers occupied at time t ∈ Z + [36], which can be written as where A t,i denotes the number of arrivals at time t − i that are still in service at time t. Note that A t,i counts how many active customers have age i. For any fixed t, (A t,i ) i≥1 , form a collection of independent and identically distributed Poisson random variables with parameter λP(S ≥ i), where λ is the average arrival rate. It is easy to calculate the expectation and variance if X t : The discrete-time process (X t ) t≥0 is stationary and time reversible. Its covariance function is Additionally, by (18), knowledge of the covariance function uniquely gives the probability mass function of the service time S because By (19), the autocovariance is non-negative and convex. Conversely, if a real-valued sequence γ k is non-negative, decreasing, and integer-convex, then it is a valid autocovariance sequence for some M/G/∞ discrete-time process [34]. If the latter assumptions hold, then lim k→∞ γ k = 0, and the probability mass function of S is exactly (19).
(X t ) t≥0 is stationary and ergodic under two assumptions: A1 At start t = 0, the number A 0,0 of customers in the system follows a Poisson probability mass function with expected value λE(S). A2 These A 0,0 customers have a service timeŜ given by a distribution which is that of the residual (or excess) life of S.
We have, as a consequence of this initial state, that (i) X t is Poisson for all t, with an expectation equal to λE(S), and (ii) the covariance function is γ k = γ 0 P( S > k) ∀k ≥ 0. Therefore, (X t ) t≥0 is LRD when S has infinite variance, for instance, whenever S has some specific discrete heavy-tailed distributions. The latter means that P(S > k) ∼ ck −q as k → ∞.

The gfGn Process
Denote by Ψ the service time in a M/G/∞ system, whose occupancy process has the covariance structure of the gfGn process. We know that at any time t, the distribution of Ψ is Poisson and converges weakly to a normal distribution when the mean goes to ∞. Specifically, combining (18) and (19), the distribution function is and the expected value is The distribution function of the residual life is needed in order to initialize the process in a steady state. It can be obtained from (18) and (20) In Figure 3, we plot the analytical distribution function of the random variable Ψ for different values of the parameters.

The Generalized Cauchy Process
Now, let Ξ be the service time in a M/G/∞ system with an occupancy process that possesses the covariance structure of the generalized Cauchy process. By (11) and (18), the cumulative distribution function of Ξ is given by and the expected value E(Ξ) is For the excess life, the distribution function is, in this case, In Figure 4, we can see the analytical distribution function of the random variable Ξ for different values of the parameters.

Accuracy
Since both distributions (21) and (24) have sub-exponential decay, the method that we developed in [37] can be used. That technique relies on the elementary properties of Poisson processes (decomposition and lack of memory) so that the M/G/∞ state process can be simulated very fast and sequentially. Moreover, the technique is flexible, i.e., valid for any distribution of the service time with sub-exponential decay. Specifically, [37] presents a fast tabular inversion method for generating the samples of S i , the service times in the queuing systems. When this tabular inversion is combined with the geometric (memoryless) interarrival times, the system state (16) can be calculated with only a few elementary operations, without any function evaluations.
The tabular method used for inversion in [37] produces some loss of numerical precision in the tail of some heavy-tailed distributions. Ref. [54] describes an improvement of the original method to solve this issue and also an extension to handle random variables with long tails in both sides. In this way, the Poisson arrival process in the case of high mean values can be also be generated using a better tabular method. Figures 5 and 6 show the estimated autocorrelation function of several traces generated with the proposed method for different values of the parameters of the generalized fGn process and the generalized Cauchy process, respectively. We can see that the estimations match well the analytical functions. be simulated very fast and sequentially. Moreover, the technique is flexible, i.e., valid for ½ any distribution of the service time with sub-exponential decay. Specifically, [37] presents ½ a fast tabular inversion method for generating the samples of S i , the service times in the ½ queueing systems. When this tabular inversion is combined with the geometric (memory-½ less) interarrival times, the system state (16) can be calculated with only a few elementary ½ operations, without any function evaluations.

½ ¼
The tabular method used for inversion in [37] produces some loss of numerical pre-½ ½ cision in the tail of some heavy-tailed distributions. [54] describes an improvement of the ½ ¾ original method to solve this issue, and also an extension to handle random variables with ½ ¿ long tails in both sides. In this way, the Poisson arrival process in the case of high mean ½ values can be also be generated using a better tabular method. ½ Figures 5 and 6 show the estimated autocorrelation function of several traces gener-½ ated with the proposed method, for different values of the parameters of the generalized ½ fGn process and the generalized Cauchy process, respectively. We can see that the estima-½ tions match well the analytical functions.

¾¼¼
In this section, we explain a method to estimate the parameters of the previous pro-¾¼½ cesses, gFGN and gCauchy, that give rise to a better adjustment of the spectral density, ¾¼¾ and therefore of the correlation structure, of empirical traffic. The method is based on the ¾¼¿ Whittle estimator [41,42], that we describe briefly in the next section. be simulated very fast and sequentially. Moreover, the technique is flexible, i.e., valid for ½ any distribution of the service time with sub-exponential decay. Specifically, [37] presents ½ a fast tabular inversion method for generating the samples of S i , the service times in the ½ queueing systems. When this tabular inversion is combined with the geometric (memory-½ less) interarrival times, the system state (16) can be calculated with only a few elementary ½ operations, without any function evaluations.

½ ¼
The tabular method used for inversion in [37] produces some loss of numerical pre-½ ½ cision in the tail of some heavy-tailed distributions. [54] describes an improvement of the ½ ¾ original method to solve this issue, and also an extension to handle random variables with ½ ¿ long tails in both sides. In this way, the Poisson arrival process in the case of high mean ½ values can be also be generated using a better tabular method. ½ Figures 5 and 6 show the estimated autocorrelation function of several traces gener-½ ated with the proposed method, for different values of the parameters of the generalized ½ fGn process and the generalized Cauchy process, respectively. We can see that the estima-½ tions match well the analytical functions.

¾¼¼
In this section, we explain a method to estimate the parameters of the previous pro-¾¼½ cesses, gFGN and gCauchy, that give rise to a better adjustment of the spectral density, ¾¼¾ and therefore of the correlation structure, of empirical traffic. The method is based on the ¾¼¿ Whittle estimator [41,42], that we describe briefly in the next section. Figure 6. Estimated autocorrelation function. Generalized Cauchy process with α = 0.5 (left) and β = 0.5 (right).

Modeling Eempirical Traces
In this section, we explain a method to estimate the parameters of the previous processes, gFGN and gCauchy, that give rise to a better adjustment of the spectral density and therefore of the correlation structure of the empirical traffic. The method is based on the Whittle estimator [41,42] that we describe briefly in the next section.

Whittle's Estimator
Let f θ (λ) be the spectral density function of a zero-mean stationary Gaussian stochastic process X = (X n ) n≥0 , where θ = (θ 1 , . . . , θ M ) is a vector of unknown parameters to be estimated from observations. Let be the empirical periodogram obtained using N samples of the process X, where · is the Euclidean norm. The Whittle estimator [41] is the vector θ = θ 1 , . . . , θ M that minimizes the statistic Whittle's estimator (28) approximates the Gaussian log-likelihood by means of the periodogram, using the latter in place of the unknown spectral density. Another informationtheoretic interpretation is that Q X (θ) is the divergence between the periodogram I X (·) and f θ (·).
Whittle's estimator converges in probability to the true value θ o . It holds that for any > 0, so θ is a weakly consistent estimator. A second property is that it is also asymptotically normal; ζ is a zero-mean Gaussian vector with a matrix of covariances known. This asymptotic normality allows the computation of confidence intervals for the estimated values. It is also possible to choose a special scale parameter such that f θ (ω) = θ 1 f * ν (ω) and the second term in (28) is zero. In this case, σ 2 is the optimal one-step-ahead prediction error that is equal to the variance of the innovations in the AR(∞) representation of the process, and Whittle's function (28) simplifies to The approximate Whittle estimator is usually evaluated numerically via the integral quadrature, replacing the integral by a sum over a discrete set of Fourier frequencies ω k = kπ n ; k = −n, . . . , −1, 0, 1, . . . , n. In addition, it can be shown that the mean squared error is σ 2 = Q X (ν), where ν is the minimizer of (31).
In the past, we proposed to use σ 2 as a measure of the suitability of a model because smaller values mean more accurate adjustment to the actual covariance function of the sample and suggest a close match between the sample and the model. Moreover, we used this value to compare the accuracy or convenience of different models. Specifically, in [37], we presented a method for model selection based on fitting the spectral density and, therefore, the empirical correlation function, and in [55], we tested the method numerically for several classes of stochastic processes, namely Gaussian and non-Gaussian processes with LRD, non-stationary processes and non-linear heteroscedastic models.

Examples
In this section, we show several examples that illustrate the suitability of these processes to model different kinds of network traffic. In these examples, we consider empirical traces that usually are used as examples of correlated traffic.
The capture of the trace of the first example was performed at the beginning of the last decade of the previous century [56]. This trace contains 2000 samples of Ethernet data lengths. The sample mean and variance are, approximately, µ = 380 and σ 2 = 180,100. In this case, both gFGN and gCauchy processes are able to approximate this trace quite well. The estimations of the parameters of each process, obtained via the Whittle estimator, are as follows: 1.
gGauchy process: β = 0.503 and α = 0.321. Figure 7 shows the correlation structure of this empirical trace versus the analytical autocorrelation of each process, with the estimated parameters.    The capture of the trace of the second example was performed in the middle of the first decade of this century [57]. This trace contains 2000 samples of TCP data lengths. The sample mean and variance are, approximately, µ = 136 and σ 2 = 46,750. In this case, only the gFGN process approximates this trace well. The estimators of its parameters, obtained via the Whittle estimator, are H = 0.661 and a = 0.961. Figure 8 shows the correlation structure of this empirical trace versus the analytical autocorrelation of the gFGN process, with the estimated parameters.  As a third example, we consider one compressed VBR video trace generated at the beginning of the second decade of this century [58]. This trace contains 2000 samples of GoPs sizes of MPEG-4 encoded video. The sample mean and variance are, approximately, µ = 34,200 and σ 2 = 5.421 × 10 8 . In this case, only the gCauchy process approximates this trace well. The estimators of its parameters, obtained via the Whittle estimator, are β = 0.221 and α = 1.962. Figure 9 shows the correlation structure of this empirical trace versus the analytical autocorrelation of the gCauchy process, with the estimated parameters. = 0.221 and α = 1.962. Figure 9 shows the correlation structure of this empirical ersus the analytical autocorrelation of the gCauchy process, with the estimated par rs.

. Discussion
In this paper, we have proposed accurate and efficient generators of the genera ractional Gaussian noise process and the generalized Cauchy process, two gener ons of previous processes that have been proposed and analysed in the last decad ood candidates for modeling different types of empirical traces. From the correl unctions, we have obtained the distribution functions of the random variables of th ice time of the M/G/∞ process, and the distribution functions of the residual life, ne order to initialize the processes in steady state. Some proofs presented show the acy of the proposed generators. Finally, we have explained how we can use the Wh stimator in order to choose the parameters of each model that give rise to a bette stment of empirical traces and for comparison of the accuracy and suitability of odel.
uthor Contributions: Both authors contributed equally to this work.

Discussion
In this paper, we proposed accurate and efficient generators of the generalized fractional Gaussian noise process and the generalized Cauchy process, two generalizations of previous processes that were proposed and analyzed in the last decade as good candidates for modeling different types of empirical traces. From the correlation functions, we obtained the distribution functions of the random variables of the service time of the M/G/∞ process and the distribution functions of the residual life needed in order to initialize the processes in steady state. Some proofs presented show the accuracy of the proposed generators. Finally, we explained how we can use the Whittle estimator in order to choose the parameters of each model that give rise to a better adjustment of the empirical traces and for comparison of the accuracy and suitability of each model.