Introduction and related works

The problem of telecommunication traffic prediction is important for telecommunications. For example, it is important for the prevention of network congestion and for the maximization of the network utilization rate [1]; it is significant for understanding future market dynamics and reducing the decision risks [2]. The telecommunication traffic prediction is also important for cyber security [3] because the comparison of the actual traffic with the predicted one may help to detect cyber-attacks.

There are a lot of different approaches to traffic prediction. For example, the following ones can be indicated: Auto Regressive Integrated Moving Average (ARIMA), Markov Modulated Poisson Process models (MMPP), Kalman filtering, Seasonal ARIMA (SA), a neural network approach (including deep neural networks [4]), wavelet transforms [1], the least-squares support vector machine (LSSVM), gray models [2], Holt-Winters models [3]. Of course, rather complicated approaches should be used for non-stationary randomly fluctuating traffic prediction. But if the traffic is stationary and rather smooth, sophisticated approaches may not be needed. For example, in [2] some methods are presented for a description of rather simple cases. In [2] it is stressed that in stationary cases the ARMA approach may be used too, and in the case of a smooth monotone process the gray model may be applied.

As is known [5], such a simple filter as the Kolmogorov-Wiener one may be used for the prediction of stationary random processes. However, as far as we know, such an approach is not sufficiently developed in the literature for traffic prediction even for rather simple cases. The Kolmogorov-Wiener filter is widely used for signal extraction in different fields of knowledge [6]. It is widely used in econometric analyses [7,8] and in image restoration [9]. The theoretical fundamentals of the Kolmogorov-Wiener filter for continuous telecommunication traffic prediction are developed in our recent paper [10]. The paper [10] is dedicated to the solution of the Wiener-Hopf integral equation in the unknown filter weight function for two telecommunication traffic models: the power-law structure function model and the model of fractional Gaussian noise; the solutions based on the truncated polynomial expansion method and the truncated trigonometric Fourier series method are obtained.

However, the possibility of using the Kolmogorov-Wiener filter for practical traffic prediction is still under question. The aim of this work is to show that the Kolmogorov-Wiener filter may be applicable to traffic prediction if the traffic is stationary and smooth enough. As is known [11,12], the telecommunication traffic in systems with data packet transfer is considered to be a self-similar heavy-tail random process. So, if we show that the Kolmogorov-Wiener filter is applicable to the prediction of simulated data of a stationary random self-similar heavy-tail process, then we will be able to conclude that it may be applied to practical telecommunication traffic prediction. In this paper we restrict ourselves to the investigation of a discrete process and a discrete filter. The corresponding simulated data may be generated via the symmetric moving average approach [13], the generated process is in fact similar to the fractional Gaussian noise process, which may describe telecommunication traffic, see [14].

The paper is organized as follows. In Sec. 1 the introduction and the literature review are given. In Sec. 2 the discrete Kolmogorov-Wiener filter and the symmetric moving average approach for obtaining simulated stationary heavy-tail data are described. In Sec. 3 heavy-tail simulated data are obtained. In Sec. 4 the prediction results are described, and in Sec. 5 conclusions are made.

Description of the discrete Kolmogorov-Wiener filter and of the method of generation of heavy-tail simulated data

Let the filter input 𝑥 𝑡 be a stationary random process which is the sum of the signal 𝑠 𝑡 and the noise 𝑛 𝑡 :

𝑥′ 𝑡 = 𝑠 𝑡 + 𝑛 𝑡 .

(1) The Kolmogorov-Wiener filter output 𝑦 𝑡 should be «the closest» to the value 𝑠 𝑡+𝑧 where 𝑧 is the number of points for which the prediction is made, so we have the following requirement:

〈(𝑦 𝑡 − 𝑠 𝑡+𝑧 ) 2 〉 → min.

(2) The correlation function 𝑅 𝑥′ (𝑡) of the filter input 𝑥′ 𝑡 and the cross-correlation function of the processes 𝑠 𝑡 and 𝑥′ 𝑡 𝑅 𝑠𝑥′ (𝑡) are considered to be given. The Kolmogorov-Wiener filter is considered to be a linear one, so the filter output is expressed in terms of the filter input as follows:

𝑦 𝑡 = ∑ ℎ 𝑖 𝑥′ 𝑡−𝑖 𝑇 𝑖=0 (3)

where ℎ 𝑖 are the unknown filter weight coefficients and the input data are given for 𝑡 = 0,1,2, . . , 𝑇. The coefficients ℎ 𝑖 should minimize expression (2). The term 〈𝑠 𝑡+𝑧 2 〉 is a constant and does not depend on the weight coefficients ℎ 𝑖 , so (2) can be rewritten as The function 𝑓(ℎ 0 , ℎ 1 , … , ℎ 𝑇 ) is a quadratic one, and thus it has one minimum, which is described by the conditions 𝜕𝑓(ℎ 0 , ℎ 1 , … , ℎ 𝑇 ) 𝜕ℎ 𝑘 = 0; 𝑘 = 0,1,2, … , 𝑇.

〈𝑦 𝑡 2 〉 − 2〈𝑦 𝑡 𝑠 𝑡+𝑧 〉 → min,(4

These conditions with account for the evenness of the correlation function and the fact that

𝜕ℎ 𝑖 𝜕ℎ 𝑗 = 𝛿 𝑖𝑗 = { 1, 𝑖 = 𝑗 0, 𝑖 ≠ 𝑗 (10) lead to ∑ ℎ 𝑖 𝑅 𝑥′ (𝑖 − 𝑘) 𝑇 𝑖=0 = 𝑅 𝑠𝑥′ (𝑘 + 𝑧); 𝑘 = 0,1,2, … , 𝑇,(11)

which is a set of linear equations in the unknown coefficients ℎ 𝑖 . In matrix form, this set may be presented as

𝑅 𝑥′ • ℎ = 𝑅 𝑠𝑥′(12)

where

𝑅 𝑥′ = ( 𝑅 𝑥′ (0) 𝑅 𝑥′ (1) 𝑅 𝑥′ (2) ⋮ 𝑅 𝑥′ (𝑇) 𝑅 𝑥′ (1) 𝑅 𝑥′ (0) 𝑅 𝑥′ (1) ⋮ 𝑅 𝑥′ (𝑇 − 1) 𝑅 𝑥′ (2) 𝑅 𝑥′ (1) 𝑅 𝑥′ (0) ⋮ 𝑅 𝑥′ (𝑇 − 2) ⋯ ⋯ ⋯ ⋱ ⋯ 𝑅 𝑥′ (𝑇) 𝑅 𝑥′ (𝑇 − 1) 𝑅 𝑥′ (𝑇 − 2) ⋮ 𝑅 𝑥′ (0) )(13)

is the correlation matrix [5], ℎ is the vector column of the unknown weight coefficients, and 𝑅 𝑠𝑥 is the vector column of the free terms:

ℎ = ( ℎ 0 ℎ 1 ℎ 2 ⋮ ℎ 𝑇) , 𝑅 𝑠𝑥′ = ( 𝑅 𝑠𝑥′ (𝑧) 𝑅 𝑠𝑥′ (𝑧 + 1) 𝑅 𝑠𝑥′ (𝑧 + 2) ⋮ 𝑅 𝑠𝑥′ (𝑧 + 𝑇)) . (14)

So, the vector column ℎ may be found as ℎ = 𝑅 𝑥′ −1 • 𝑅 𝑠𝑥′ .

(15) Then the filter output may be obtained by formula (3).

It should be noticed that all the above-mentioned calculations are described in [6]. The Kolmogorov-Wiener filter may be used both for the extraction of a signal form the sum of a signal and a noise and for the signal prediction. In the case where the input signal is non-noisy, the Kolmogorov-Wiener filter may be used for the prediction of the stationary process given at the filter input. In the non-noisy case, the filter weight coefficients are given by formula (15) with account for the fact that 𝑅 𝑠𝑥′ = (𝑅 𝑥′ (𝑧) 𝑅 𝑥′ (𝑧 + 1) 𝑅 𝑥′ (𝑧 + 2) … 𝑅 𝑥′ (𝑧 + 𝑇)) 𝑇 .

(16) Now let us describe the method of the generation of heavy-tail simulated data, which is used in the paper. In this paper we use the symmetric moving average approach, which is described in detail in [13]. Such an approach was chosen because of its simplicity.

Let 𝑉 𝑡 be a stationary white noise process with an average value equal to zero and a variance equal to 1. Then a heavy-tail process 𝑋 𝑖 similar to the fractional Gaussian noise may be generated as follows [13]:

𝑋 𝑖 = ∑ 𝑎 |𝑗| 𝑉 𝑖+𝑗 𝑞 𝑗=−𝑞 = 𝑎 𝑞 𝑉 𝑖−𝑞 + 𝑎 𝑞−1 𝑉 𝑖−𝑞+1 + ⋯ + 𝑎 𝑞 𝑉 𝑖+𝑞 ,(17)

theoretically, 𝑞 should be infinite; in practical calculation it may be a rather large, but finite number; and the coefficients 𝑎 𝑗 are as follows:

𝑎 0 = √(2 − 2𝐻)𝛾 0 1.5 − 𝐻(18)

and

𝑎 𝑗 = 𝑎 0 2 ((𝑗 + 1) 𝐻+0.5 + (𝑗 − 1) 𝐻+0.5 − 2𝑗 𝐻+0.5 ),(19)

here, 𝛾 0 is the variance and 𝐻 is the Hurst exponent of the process 𝑋 𝑖 . The number 𝑞 may be very large, it is estimated as follows [13]:

𝑞 ≥ max (𝑚, ( 𝐻 2 − 0.25 2𝛽 ) 1 1.5−𝐻 ) (20)

where 𝑚 is the number of correlation function points of the process 𝑋 𝑖 which should be obtained and a small number 𝛽 is in fact the given accuracy of the coefficient 𝑎 𝑗 in (17); the values 𝑎 𝑗>𝑞 should be less than 𝛽𝑎 0 . The accuracy of this method depends on 𝑞, and the method is not exact even in the case where 𝑞 → ∞. However, for a rather large 𝑞 the method may lead to good practical results [13].

3. The generation of non-smooth and smooth heavy-tail simulated data 10 6 points of the white noise process 𝑉 𝑡 with an average value equal to 0 and a variance equal to 1 are generated on the basis of the generator built in the Wolfram Mathematica package. The following parameters were chosen: 𝑚 = 10 5 , 𝛽 = 10 −4 , 𝐻 = 0.8, 𝛾 0 = 1.

(21) The corresponding number 𝑞 = 3 • 10 5 is chosen. In fact, the inequality (21) holds even for 𝑞 = 10 5 , the value 𝑞 = 3 • 10 5 was chosen for a higher accuracy. On the basis of the idea (17)-(19), 10 5 points of the process 𝑋 𝑖 were generated as follows:

𝑋 𝑖 = ∑ 𝑎 |𝑗| 𝑉 𝑖+𝑗+𝑞 𝑞 𝑗=−𝑞 , (22)

in fact, the quantities 𝑉 𝑖+𝑗+𝑞 and 𝑉 𝑖+𝑗 are independent because 𝑉 𝑡 is the white noise, no matter whether formula (17) or formula (22) is used; formula (22) is chosen in order to avoid indices beyond the array 𝑉 𝑖 bounds. The coefficients 𝑎 𝑗 are calculated on the basis of (19).

The average value of 𝑋 𝑖 is close to zero. We have to construct simulated data that may describe telecommunication traffic, which is obviously non-negative. So we build the array 𝑥 𝑖 as follows:

𝑥 𝑖 = 𝑋 𝑖 + |min(𝑋)| + 10 −3 , (23) a small summand 10 -3 is added in order to avoid obtaining an infinite value of the prediction mean average percentage error (MAPE). The process 𝑥 𝑖 is a non-negative random stationary heavy-tail process; its graph is given in Fig. 1.

Let us make sure that the generated process 𝑥 𝑖 is a heavy-tail one. Let us describe the corresponding centralized process 𝑥𝑐 𝑖 :

𝑥𝑐 𝑖 = 𝑥 𝑖 − 〈𝑥〉 (24) where the average value 〈𝑥〉 is

〈𝑥〉 = 1 10 5 ∑ 𝑥 𝑖 10 5 𝑗=1 ,(25)

here we take into account the fact that the number of points of the generated array 𝑥 𝑖 is equal to 10 5 . The correlation function of the process 𝑥𝑐 𝑖 is built as follows:

𝑅 𝑥 (𝜏) = 〈𝑥𝑐 𝑖 • 𝑥𝑐 𝑖+𝜏 〉 = 1 10 5 − 𝜏 ∑ (𝑥𝑐 𝑖 • 𝑥𝑐 𝑖+𝜏 ) 10 5 −𝜏 𝑖=1 . (26)

The corresponding correlation function and its least-square fit are given in Fig. 2. The corresponding least-square fit is sought as 𝑅 fit (𝑡) = 𝑎 • 𝑡 𝑏 , (27) The following numerical coefficients were obtained: 𝑎 = 0.39, 𝑏 = −0.44, (28) here, the coefficients are rounded off to two significant digits. So, 𝑅 fit (𝑡) = 0.39 • 𝑡 −0.44 , (29) and on the basis of formula (29) and Fig. 2 one can conclude that the correlation function exhibits a power law decay rather than an exponential one. So, indeed, the generated process is a heavy-tail one.

It should also be noticed that according to [13] the following property should be valid for ≥ 1 : 𝑅 𝑥 (𝑡)~𝑡 2𝐻−2 , (30) so, according to the least-square fit 2𝐻 − 2 = −0.44, (31) which leads to 𝐻 = 0.78, (32) which is very close to the value 0.8, see (21). The variance of the process is equal to 𝑅 𝑥 (0) = 0.93, (33) which is rather close to the value 𝛾 0 = 1, see (21). So one can conclude that the generated process is close to the fractional Gaussian noise with given variance and Hurst exponent.

The generated process is non-smooth, i.e. it is really highly fluctuating, so it is rather difficult to predict it. So it is reasonable to investigate smooth heavy-tail processes. In order to obtain smoother processes, we use a very simple smoothing algorithm [15]:

𝑋 ̃𝑖 = 1 2𝑙 + 1 ∑ 𝑋 𝑖+𝑗 𝑙 𝑗=−𝑙 (34)

where 𝑋 ̃𝑖 are the values of a smooth process, expression (34) is valid for every point except for the first 𝑙 and the last 𝑙 ones. The first 𝑙 and the last 𝑙 points of the process 𝑋 ̃𝑖 may be obtained as the corresponding linear least-square fit of the first 𝑙 and the last 𝑙 points of the process 𝑋 𝑖 , respectively. The corresponding non-negative process may be expressed similarly to ( 23 .

(37)

The simulated data for the process 𝑥 ̃𝑖 for 𝑙 = 3 are given in Fig. 3. It should be stressed that the obtained smooth process 𝑥 ̃𝑖 is also a heavy-tail one. Let us consider the corresponding correlation function:

𝑅 𝑥 ̃(𝜏) = 〈𝑥𝑐 ̃𝑖 • 𝑥𝑐 ̃𝑖+𝜏 〉 = 1 10 5 − 𝜏 ∑ (𝑥𝑐 ̃𝑖 • 𝑥𝑐 ̃𝑖+𝜏 ) 10 5 −𝜏 𝑖=1 . (38)

For example, for 𝑙 = 3 the following correlation function and its fit are obtained, see Fig. 4. The least-square fit is sought in the form (27) , the following numerical coefficients were obtained: 𝑎 = 0.43, 𝑏 = −0.46, (39) here, the coefficients are rounded off to two significant digits. So, 𝑅 fit (𝑡) = 0.43 • 𝑡 −0.46 . (40) As can be seen form Fig. 4, the correlation function of a smooth process is also well described by a power-law function, the obtained smooth process 𝑥 ̃𝑖 is also a heavy-tail one, and, in fact, this process may also be roughly considered as fractional Gaussian noise.

Prediction on the basis of the Kolmogorov-Wiener filter

The prediction for non-smooth data is built as follows. In fact, the prediction for the centralized process is used. The filter weight coefficients are built on the basis of ( 13)-( 16); the corresponding correlation function is taken in the form (26).

First of all, the points 𝑥𝑐 1 , 𝑥𝑐 2 ,…, 𝑥𝑐 𝑇+1 of the simulated process 𝑥𝑐 are taken as the filter input, and the points 𝑥𝑐 𝑇+2 , 𝑥𝑐 𝑇+3 ,…, 𝑥𝑐 𝑇+𝑧+1 are predicted. Then the points 𝑥𝑐 2 , 𝑥𝑐 3 , … , 𝑥𝑐 𝑇+2 are taken from the simulated data, and the points 𝑥𝑐 𝑇+3 , 𝑥𝑐 𝑇+4 ,…, 𝑥𝑐 𝑇+𝑧+2 are predicted, and so on.

At the 𝑖 th iteration of the algorithm the predcition is calculated as follows. The filter input data are

𝑥′ 0 = 𝑥𝑐 𝑖 , 𝑥′ 1 = 𝑥𝑐 𝑖+1 , … , 𝑥′ 𝑇 = 𝑥𝑐 𝑖+𝑇 ,(41) so

𝑥′ 𝑗 = 𝑥𝑐 𝑖+𝑗 .

(42) The filter output 𝑦 𝑡 is the predicted value for 𝑥′ 𝑡+𝑧 (the non-noisy case is investigated). According to (3) we have

𝑦 𝑡 = ∑ ℎ 𝑘 𝑥′ 𝑡−𝑘 , 𝑡 𝑘=0(43)

here, the upper bound of summation is changed in order to avoid obtaining indices beyond the array of the input data. Such a change of the bound does not lead to a significant error for the prediction under consideration. On the basis of (41)-(43) one can conclude that

𝑥𝑐 ̂𝑖+𝑡+𝑧 = ∑ ℎ 𝑘 𝑥𝑐 𝑡+𝑖−𝑘 𝑡 𝑘=0 (44)

where 𝑥𝑐 ̂𝑖+𝑡+𝑧 is the predicted value of 𝑥𝑐 𝑖+𝑡+𝑧 . Obviously, the prediction is made only for the values 𝑖 + 𝑡 + 𝑧 = 𝑇 + 1 + 𝑖, 𝑇 + 2 + 𝑖, … , 𝑇 + 𝑧 + 𝑖. We should also remember that we should make the prediction for the non-negative simulated data. So, the predicted non-negative data may be expressed as

𝑥 ̂𝑖+𝑡+𝑧 = 𝑥𝑐 ̂𝑖+𝑡+𝑧 + 〈𝑥〉 = 〈𝑥〉 + ∑ ℎ 𝑘 𝑥𝑐 𝑡+𝑖−𝑘 , 𝑡 = 𝑇 + 1 − 𝑧, 𝑇 ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ 𝑡 𝑘=0 .(45)

The MAPE and MAE errors for the corresponding prediction are calculated as

MAPE = 1 𝑧 ∑ | 𝑥 ̂𝑖+𝑡+𝑧 − 𝑥 𝑖+𝑡+𝑧 𝑥 𝑖+𝑡+𝑧 | 𝑇 𝑡=𝑇+𝑧−1 • 100%(46)

and

MAE = 1 𝑧 ∑ |𝑥 ̂𝑖+𝑡+𝑧 − 𝑥 𝑖+𝑡+𝑧 | 𝑇 𝑡=𝑇+𝑧−1 .(47)

The corresponding prediction errors are calculated at each iteration. Let us tell a few words why the above-mentioned change of the upper bound of summation has no significant effect on the result. In order to make the prediction for 𝑥𝑐 ̂𝑇+1+𝑖 , one should calculate the sum of 𝑇 + 2 − 𝑧 summands, in order to make the prediction for 𝑥𝑐 ̂𝑇+2+𝑖 one should calculate the sum of 𝑇 + 3 − 𝑧 summands, and so on. We obviously deal with the case where 𝑇 ≫ 𝑧, so the value 𝑇 + 1 − 𝑧 is rather close to 𝑇 + 1, so the above-mentioned change of the upper bound is not significant for the calculations.

Similarly, the prediction for the smooth heavy-tail process is made as follows. At the 𝑖 th iteration of the algorithm the prediction is calculated as follows: The following results are obtained. The MAPE and MAE histograms in the case of the non-smooth process are shown in Fig. 5. The y-axes of the histograms indicate the number of MAPE and MAE values that belong to the corresponding intervals. For the non-smooth process the average MAPE is 24.7%, and the average MAE is 0.70 (the average value of the process is 〈𝑥〉 = 3.88). It should also be stressed that for some points the MAPE are more than 100%. So one can conclude that the prediction accuracy is not high in the case of the non-smooth process. So, if the process is a highly fluctuating one, then the prediction based on the Kolmogorov-Wiener filter may not lead to good results.

𝑥

But if the process is rather smooth, the prediction results are much better. The corresponding results are given in Table 1. In Table 1 𝑙 is the parameter used in (34), i.e. 2𝑙 + 1 is the number of smoothing points. As can be seen, the smother the process is, the better the prediction results are, and the prediction accuracy increases with 𝑙. The corresponding histograms for 𝑙 = 3 are given in Fig. 6. The predictions for 𝑙 ≥ 6 have an average MAPE value less than 3%.

Figure 6:

The MAPE and MAE histograms for the prediction of a non-smooth heavy-tail process For example, for 𝑙 = 3 the average MAPE is less than 5%. As can be seen from the corresponding histogram, the MAPE for the overwhelming majority of points is less than 10%. For some very rare points the MAPE may be rather high (up to 40%), but in our opinion this may be explained as follows. As can be seen from Fig. 3, the values for some points of the process 𝑥 ̃ are rather close to zero, and the MAPE may not be an adequate characteristic for the prediction of points close to zero. So, one can conclude that the Kolmogorov-Wiener filter may give good results for the prediction of a stationary heavy-tail random process if the process is smooth enough.

Conclusions and plans for the future

The use of the Kolmogorov-Wiener filter for the prediction of stationary random heavy-tail processes is considered. The attention is paid to the discrete case. The problem under consideration may be connected with the telecommunication traffic prediction, which is important, for example, for cyber security, see [3]. There are many rather sophisticated approaches to telecommunication traffic prediction [1]. For rather simple cases (stationary or smooth traffic) the ARMA or gray model approaches may be used [2]. The traffic in telecommunication systems with data packet transfer is considered to be a self-similar heavy-tail process, see [11]. Such a simple filter as the Kolmogorov-Wiener one may be used in the prediction of stationary random processes [6]. However, as far as we know, the corresponding approach for traffic prediction is not sufficiently developed in the literature.

In this paper we generate data for a stationary heavy-tail process on the basis of the symmetric moving average approach [13]. The corresponding non-smooth and smooth data are generated. The prediction for 1 point forward on the basis of the previous 101 points is investigated. It is shown that the Kolmogorov-Wiener filter is not good for non-smooth processes, but may give a good prediction for a stationary random heavy-tail process if the process is rather smooth. So, if the traffic is stationary and rather smooth, the Kolmogorov-Wiener filter may be used for its prediction. The advantage of the corresponding approach is the simplicity of the method in contrast with, for example, neural networks or ARIMA models.

The plans for the future are as follows. In this paper only the values T = 100 and z = 1 are investigated. So the prediction investigation for a wider range of parameters may be a plan for the future. In our recent paper [10] the theoretical approach to the Kolmogorov-Wiener filter construction in the continuous case is considered. In this paper we generated a large number of data points, which may allow one to try to investigate the continuous case, so the investigation of the applicability of the method [10] may be another plan for the future. This paper is based on the generation of simulated data, so the investigation of real experimental traffic data may be another plan for the future. It should also be stressed that the use of the Kolmogorov-Wiener filter for the prediction of stationary