Novel anomaly detection approach for telecommunication network proactive performance monitoring

Abstract

The mode of telecommunications network management is changing from “network oriented” to “subscriber oriented”. Aimed at enhancing subscribers’ feeling, proactive performance monitoring (PPM) can enable a fast fault correction by detecting anomalies designating performance degradation. In this paper, a novel anomaly detection approach is the proposed time series prediction and the associated confidence interval based on multiplicative autoregressive integrated moving average (ARIMA). Furthermore, under the assumption that the training residual is a white noise process following a normal distribution, the associated confidence interval of prediction can be figured out under any given confidence degree 1-α by constructing random variables satisfying t distribution. Experimental results verify the method’s effectiveness.

Keywords

proactive performance monitoring (PPM) / anomaly detection / time series prediction / autoregressive integrated moving average (ARIMA) / white noise / confidence interval

Introduction

The mode of telecommunication network management is changing from “network oriented” to “subscriber oriented”. Accordingly, the mode of performance management of network, an important aspect of network management, is transforming from passiveness to initiativeness. The proactive performance monitoring (PPM) is simply one of the initiative management ways.

PPM [1,2] is concerned with performance analysis and fault detection, which are capable of detecting “soft” network and service faults automatically and adaptively in the midst of networks’ performance fluctuations and evolutions. By detecting anomalies designating performance degradation, which are the symptoms of network faults and preludes to services failure, PPM can enable a fast fault containment and correction, through which serious network failures can be avoided and duration time can be shortened. “Anomaly” is defined as “statistically unusual”. It can be detected by using two terms: baseline and thresholds (upper and lower thresholds). Thresholds fluctuate around predictive baseline. If the actual value of key performance indicator (KPI) exceeds thresholds, performance warning would be triggered to show the deteriorating performance quality.

For anomaly detection, the prediction of the baseline and thresholds of KPI is a crucial point. If the threshold is set over strictly, a fake warning may be launched. However, if the threshold is too high or too low, the mechanism will make little sense. There were some researches on the derivations of baseline and thresholds [2,3]. In Ref. [2], the threshold that is called tolerance limits or envelope is based on standard deviation of training set, while in Ref. [3], the baseline is based on a simplistic algorithm, and they are all too coarse. In this paper, time series prediction is deeply studied to settle baseline, and confidence interval of prediction error is used for thresholds. Under the assumption that a white noise process follows normal distribution, the associated confidence interval of prediction value can be worked out under any dedicated confidence degree 1-α by constructing random variables following t distribution.

This paper is organized as follows. In Sect. 2, the principle and model identification method of autoregressive integrated moving average (ARIMA) is introduced. In Sect. 3, the proposed time series prediction approach with the associated confidence interval calculation using multiplicative ARIMA is presented in particular. In Section 4, an experiment with the proposed approach is carried out. The conclusion is drawn in Sect. 5.

Overview of ARIMA

Let time series x₁,x₂,…,x_i,…,x_N denote the observations made at equidistant time interval {\tau }_{\mathrm{0}}+h, {\tau }_{\mathrm{0}}+\mathrm{2}h, …, {\tau }_{\mathrm{0}}+ih, …, {\tau }_{\mathrm{0}}+Nh, where x_i represents observation at time {\tau }_{\mathrm{0}}+ih if we consider {\tau }_{\mathrm{0}} as the origin and h as the unit of time. Basically, the time series prediction can be considered as a modeling issue; a model is built between input and output. Then, the model is used to predict the future values based on the previous values.

Among the algorithms deployed in time series prediction, both autoregressive moving average (ARMA) [4] and ARIMA [5] are most commonly used. ARMA can only be used for linear and stationary time series, while ARIMA can be used for nonstationary ones exhibiting homogeneity. When the time series show heterogeneous nonstationarity, other algorithms such as artificial neural network (ANN) [6] and support vector machine (SVM) will be adopted [7-9].

In this paper, we use multiplicative ARIMA, ARIMA with seasonal part, to model KPI time series with homogeneous nonstationarity.

Since ARIMA is a model evolving from ARMA, which can be regarded as an integration of autoregressive (AR) model and moving average (MA) model, AR, MA, and ARMA models are introduced first.

AR process

Autoregressive model has the form as follows:

(1)x_t={\displaystyle \sum _{i=\mathrm{1}}^p{\varphi }_i{x}_{t-i}+{ϵ}_t},

where {\varphi }_p\ne \mathrm{0}, and {ϵ}_t is a white noise process satisfying WN(\mathrm{0},{\sigma }^{\mathrm{2}}).

Accordingly, the process defined by Eq. (1) is called an autoregressive process of order p, or more succinctly, an AR(p) process.

MA process

Moving average model has the form as follows:

(2)Y_t={ϵ}_t-{\displaystyle \sum _{i=\mathrm{1}}^q{\theta }_i{ϵ}_{t-i}},

where {\theta }_q\ne \mathrm{0}, and {ϵ}_t is a white noise process satisfying WN(\mathrm{0},{\sigma }^{\mathrm{2}}).

The process defined by Eq. (2) is called moving average process of order q, which is abbreviated to MA(q).

Mixed ARMA process

Sometimes, a series can be modeled as autoregressive process, while some others can be modeled as moving average process. However, there are time series that cannot be modeled as pure AR or MA due to the too many parameters required. In these cases, mixed ARMA model should be used where both autoregressive and moving average terms are in the model, which can be represented as

(3)x_t={\displaystyle \sum _{i=\mathrm{1}}^p{\varphi }_i{x}_{t-i}}-{\displaystyle \sum _{j=0}^q{\theta }_j{ϵ}_{t-j}},

or

(4)\varphi (B)x_t=\theta (B){ϵ}_t.

In Eq. (4), B stands for backward shilt operator, which is defined by B{x}_t=x_{t-1},

A process defined by Eq. (3) or (4) is called mixed autoregressive moving average process of order (p, q), which is abbreviated to ARMA(p,q).

ARIMA process

AR, MA, and ARMA models described above are all appropriate for modeling stationary process. However, in reality, many time series are nonstationary. For those nonstationary time series, nevertheless exhibiting homogeneity, we can stationalize them by differencing operation. Let \nabla =\mathrm{1}-B denote the differencing operator so that ARIMA model can be written as

(5)\varphi (B){\nabla }^d{x}_t=\theta (B){ϵ}_t,

where d is an integer representing number of differencing operation, \varphi (B) is a stationary autoregressive operator, and \theta (B) is a moving average operator. Thus, we can get the conclusion that the model represented by Eq. (5) corresponds to the assumption that the dth difference of the series can be represented by a stationary ARMA process. The process satisfying Eq. (5) is called an ARIMA process.

Similarly, if the nonstationary time series exhibit periodical homogeneity with period s, they can be stationalized by seasonal differencing. The seasonal operator is of the form {\nabla }_s=\mathrm{1}-B^s. If a time series is only of seasonal nonstationary, the model can be represented as

(6){\mathrm{\Phi }}_P(B^s){\nabla }_s^D{x}_t={\mathrm{\Theta }}_Q(B^s){\omega }_t,

where

{\mathrm{\Phi }}_P(B^s)=\mathrm{1}-{\mathrm{\Phi }}_{\mathrm{1}}{B}^s-{\mathrm{\Phi }}_{\mathrm{2}}{B}^{\mathrm{2}s}-{\mathrm{\Phi }}_{\mathrm{3}}{B}^{\mathrm{3}s}-\cdots -{\mathrm{\Phi }}_P{B}^{Ps},

and

{\mathrm{\Theta }}_Q(B^s)=\mathrm{1}-{\mathrm{\Theta }}_{\mathrm{1}}{B}^s-{\mathrm{\Theta }}_{\mathrm{2}}{B}^{\mathrm{2}s}-{\mathrm{\Theta }}_{\mathrm{3}}{B}^{\mathrm{3}s}-\cdots -{\mathrm{\Theta }}_Q{B}^{Qs}.

If time series have both trendy and seasonal part, integrating Eqs. (5) and (6), we can get the following model:

(7)\varphi (B){\mathrm{\Phi }}_P(B^s){\nabla }^d{\nabla }_s^D{Y}_t=\theta (B){\mathrm{\Theta }}_Q(B^s){ϵ}_t.

The resulting multiplicative ARIMA process will be of order (p,d,q)\times {(P,D,Q)}_s.

Time series prediction based on ARIMA with associated confidence interval

For some KPIs such as traffic volume or call attempts, the baseline can be obtained according to time series prediction using ARIMA. The threshold then can be achieved based on the confidence interval under a dedicated confidence degree \mathrm{1}-\alpha by using the hypothesis that training residual satisfying white noise follows normal distribution.

Take the hourly traffic of mobile service switching center (MSC) for example. The traffics in different hours in a day are usually different. Because PPM takes hours as a unit, we predict hourly traffic using traffics at the same time in the past days. For example, if we want to predict the traffic at 9:00 tomorrow, then the history values are traffic data at 9:00 of today and yesterday and so on. Here, we use traffic values at 9:00 from 2007-7-4 to 2007-8-17, which is totally 45 consecutive items, as training data.

There are three steps to make the time series prediction using ARIMA. First, the ARIMA model should be identified. Second, the parameters should be estimated for the identified class of ARIMA models. Third, the prediction for future values should be made by using history values. In this paper, the associated confidence interval of prediction value should also be computed under given confidence degree.

The three steps above will be discussed in detail as follows.

Model identification

To get the ARIMA model fit for a time series, the order of the model (p,d,q)\times {(P,D,Q)}_s should be decided at the first stage. Thus, an appropriate subclass of models can be identified from the general seasonal multiplicative ARIMA family denoted by Eq. (7). This stage is usually called model identification. By calculating the autocorrelation function, we can justify if the series exhibit trend or if the series shows seasonal pattern with period s. Then, we can make the series stationary by using a differencing operator or seasonal differencing operator. For stationary time series, there are several model identification criteria including AIC, BIC, F-test, etc.

In an actual application, especially when the series length of available data is not very large, the orders P, Q, p, and q would rarely need to be greater than 1. Therefore, at this stage of model identification, we should select the best model from all the combinations of (p,d,q)\times {(P,D,Q)}_s where all the orders are between zero and one according to AIC.

Parameter estimation

After the appropriate model (p,d,q)\times (P,D,Q) has been identified, the parameters would be estimated. There are several estimation algorithms, such as Yule-Walker equation, minimum mean squared error (MMSE), and maximum likelihood estimation (MLE). MLE is used in this paper.

Prediction with associated confidence interval

After the model is built up, it can be used to forecast future values based on history data. There are two types of forecasting, one-step ahead forecasting, and multisteps ahead forecasting. To guarantee the precision of forecasting, one-step ahead forecasting is adopted.

Actual value at time {\tau }_{\mathrm{0}}+ih is denoted by x_i, and prediction value for the same time is denoted by {\widehat{x}}_i, the equation

x_i={\widehat{x}}_i+{ϵ}_i

holds, where {ϵ}_i is the prediction error and is white noise process following normal distribution, {ϵ}_i\sim N(\mathrm{0},{\sigma }^{\mathrm{2}}). Therefore, the actual value at the monitoring time point {\widehat{x}}_{n+\mathrm{1}} and prediction value {\widehat{x}}_{n+\mathrm{1}} satisfy {\stackrel{\square }{x}}_{n+\mathrm{1}}\sim N\left({\widehat{x}}_{n+\mathrm{1}},{\sigma }^{\mathrm{2}}\right), which is equivalent to {\displaystyle \frac{x_{n-\mathrm{1}}-{\widehat{x}}_{n+\mathrm{1}}}{\sigma }}\sim N(\mathrm{0},\mathrm{1}). Using a normal distribution table, we can compute the confidence interval of actual value 1 under a given confidence degree 1-α, and the following equation is established:

P\left\{-z_{\mathrm{\pi }/2}\le {\displaystyle \frac{x_{n+\mathrm{1}}-{\widehat{x}}_{n+\mathrm{1}}}{\sigma }}\le z_{\mathrm{\pi }/2}\right\}=\mathrm{1}-\alpha .

If the number of training items is more than 50, the training sample standard deviation

S=\sqrt{{\displaystyle \frac{\mathrm{1}}{k-\mathrm{1}}}{\displaystyle \sum _{i=\mathrm{1}}^k{\left({ϵ}_i-{\overline{ϵ}}_i\right)}^{\mathrm{2}}}}

can be regarded as equal to the population standard deviation σ. Let a confidence degree be \mathrm{1}-\alpha =\mathrm{0.95}, so that we can get the confidence interval of x_n+1 as

\left({\widehat{x}}_{n+\mathrm{1}}-\mathrm{2}\sigma ,{\widehat{x}}_{n+\mathrm{1}}+\mathrm{2}\sigma \right)=\left({\widehat{x}}_{n+\mathrm{1}}-\mathrm{2}S,{\widehat{x}}_{n+\mathrm{1}}+\mathrm{2}S\right).

If the confidence degree is \mathrm{1}-\alpha =\mathrm{0.97}, the confidence interval will be \left({\widehat{x}}_{n+\mathrm{1}}-\mathrm{3}\sigma ,{\widehat{x}}_{n+\mathrm{1}}+\mathrm{3}\sigma \right), which can be substituted with \left({\widehat{x}}_{n+\mathrm{1}}-\mathrm{3}S,{\widehat{x}}_{n+\mathrm{1}}+\mathrm{3}S\right).

However, if the number of training items is less than 50, the value of both σ and S may vary significantly so that they can not be regarded as being equivalent. In this case, we propose the following algorithm to get the confidence interval.

Since training residual \{{ϵ}_i\} satisfies normal distribution, according to the linear characteristics of operation of random variables obeying normal distribution, its mean value denoted by \overline{ϵ} also satisfies normal distribution:

\overline{ϵ}\sim N\left(\mathrm{0},{\displaystyle \frac{{\sigma }^{\mathrm{2}}}{n}}\right).

The random variable {ϵ}_{n+\mathrm{1}}-\overline{ϵ} will also satisfies the following normal distribution:

(8){ϵ}_{n+\mathrm{1}}-\overline{ϵ}\sim N\left(\mathrm{0},\left[\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{n^{\mathrm{2}}}}\right]{\sigma }^{\mathrm{2}}\right),

and Eq. (8) can be transformed to

(9)U={\displaystyle \frac{{ϵ}_{n+\mathrm{1}}-\overline{ϵ}}{\sqrt{{\displaystyle \frac{n^{\mathrm{2}}+\mathrm{1}}{n^{\mathrm{2}}}}}\sigma }}\sim N(\mathrm{0},\mathrm{1}).

Considering the fact that the standard deviation S of training residual satisfying χ² distribution, a random variable V can be derived as follows:

(10)V={\displaystyle \frac{(n-\mathrm{1})S^{\mathrm{2}}}{{\sigma }^{\mathrm{2}}}}~{\chi }^{\mathrm{2}}(n-\mathrm{1}).

Integrating Eqs. (9) and (10), another random variable named Z satisfying t distribution can be obtained by

(11)Z={\displaystyle \frac{U}{\sqrt{{\displaystyle \frac{V}{n-\mathrm{1}}}}}}=\sqrt{{\displaystyle \frac{n^{\mathrm{2}}}{n^{\mathrm{2}}+\mathrm{1}}}}{\displaystyle \frac{({ϵ}_{n+\mathrm{1}}-\overline{ϵ})}{S}}~t(n-\mathrm{1}).

In Eq. (11), only {ϵ}_{n+\mathrm{1}} is unknown, it can be transformed to

(12)P\left\{\overline{ϵ}-\sqrt{{\displaystyle \frac{n^{\mathrm{2}}+\mathrm{1}}{n^{\mathrm{2}}}}}{t}_{\mathrm{\pi }/2}(n-\mathrm{1})S\le {ϵ}_{n+\mathrm{1}}\le \overline{ϵ}+\sqrt{{\displaystyle \frac{n^{\mathrm{2}}+\mathrm{1}}{n^{\mathrm{2}}}}}{t}_{\mathrm{\pi }/2}(n-\mathrm{1})S\right\}=\mathrm{1}-\alpha .

Therefore, the confidence interval of x_n+1 under confidence degree \mathrm{1}-\alpha will be

\left[\overline{ϵ}-\sqrt{{\displaystyle \frac{n^{\mathrm{2}}+\mathrm{1}}{n^{\mathrm{2}}}}}{t}_{\mathrm{\pi }/2}(n-\mathrm{1})S+{\widehat{x}}_{n+\mathrm{1}},\overline{ϵ}+\sqrt{{\displaystyle \frac{n^{\mathrm{2}}+\mathrm{1}}{n^{\mathrm{2}}}}}{t}_{\mathrm{\pi }/2}(n-\mathrm{1})S+{\widehat{x}}_{n+\mathrm{1}}\right].

Thus, the former value is taken as the lower threshold, and the latter is the upper threshold with the confidence degree \mathrm{1}-\alpha.

Anomaly detection for traffic series using ARIMA

Now, we have an hourly traffic load series at 9:00 from 2007-7-4 to 2007-8-17 with length of 45, as shown in Table 1. From the first item, we always take 37 consecutive items as training data to predict the next item with the associated confidence interval.

We stationalize the original traffic time series by differencing operation. The autocorrelation function of the original traffic series is depicted in Fig. 1. This figure reveals that there is no trend part but a marked seasonal pattern with periodicity of 7, so a seasonal difference should be carried out and denotes the resulting series by {y_i}, where y_i={\nabla }_{\mathrm{7}}{x}_i=(\mathrm{1}-B^{\mathrm{7}})x_i. Then, the stationary of {y_i} is explored by calculating autocorrelation function. The result is shown in Fig. 2, where the seasonally differenced series { y_i } is stationary.

We model the training data using ARIMA (p,0,q)×(P,1,Q₇) model and identify the model using AIC. Then, all the combinations of (p,0,q)×(P,1,Q) where each order is 0 or 1 were tested and ensured that they are not zeros at the same time and select the combination that makes AIC the least. The model (1,0,1)×(1,1,0) is selected, which means that the original time series {x_i} (1,0,1)×(1,1,0)₇ is the best model.

We make a one-step ahead prediction as the baseline and calculate the associated confidence interval as the threshold. The one-step ahead prediction value is 2249.6, which was taken as the hourly traffic for 9:00, 2007-8-10. The sample standard deviation is S=42.03. With the confidence degree \mathrm{1}-\alpha =\mathrm{0.95}, we have

t_{\alpha /2}(n-\mathrm{1})=\mathrm{2.05},

and

\sqrt{{\displaystyle \frac{n^{\mathrm{2}}+\mathrm{1}}{n^{\mathrm{2}}}}}{t}_{\alpha /2}(n-\mathrm{1})S=\mathrm{86.212.}

With the confidence degree \mathrm{1}-\alpha =\mathrm{0.95}, the confidence interval is

[\mathrm{2249.6}+\mathrm{2.77}-\mathrm{86.212},\mathrm{2249.6}+\mathrm{2.77}+\mathrm{86.212}]=[\mathrm{2166.158},\mathrm{2338.582}].

Therefore, we reach conclusion that the actual traffic value for 2007-8-10 9:00 is 2238.39, which is among the confidence interval [2166.158, 2338.582]. This means that during this period, the network is in normal state.

Using the method above, the prediction with associated confidence interval can be made from 2007-8-10 until 2007-8-17. The result is shown in Table 2 and plotted in Fig. 3.

We come to the following conclusions based on Table 2 and Fig. 3:

1) The actual traffics at 9:00 from 2007-8-10 till 2007-8-16 are all within the tolerance limit and all the values of absolute percent errors (APEs) are smaller than 3%. Figure 3 shows that in normal state, the proposed prediction approach is of high precision.

2) However, at 9:00 on 2007-8-17, the actual traffic, 2101.47, falls outside tolerance limit [2156.72, 2364.56]. In this status, a performance alarm should be launched. This abnormal status has been confirmed by the fact that during this period, there was something wrong with IC2 module in BSU unit, which caused some call setup failure. That is why the traffic was so low then. 2007-8-17 is Friday when traffic is usually the heaviest in the week. However, Fig. 2 shows us an opposite case that the traffic on 2007-8-17 is the lightest in the week, which is obviously abnormal.

Conclusion

The prediction of baseline and threshold for KPIs is a crucial issue for proactive performance monitoring of telecommunication network. In this paper, the approach of time series prediction with the associated confidence interval using ARIMA was adopted to solve this problem. Furthermore, under the assumption that a white noise process follows normal distribution, the associated confidence interval of prediction can be computed under any given confidence degree \mathrm{1}-\alpha by constructing random variable satisfying t distribution. Experimental results show that using the proposed modeling method, a quite precise prediction and fluctuation range can be achieved.

References

[1] Hellerstein J L, Zhang F, Shahabuddin P. An approach to predictive detection for service management. In: Proceedings of the 6th IFIP/IEEE International Symposium on Integrated Network Management, 1999, Sloman M, Mazumdar S, Lupu E, Eds. New York: IEEE Publishing, 1999, 309-322

[2] Feather F, Siewiorek D, Maxion R. Fault detection in an Ethernet network using anomaly signature matching. ACM SIGCOMM Computer Communication Review, 1993, 23(4): 279-288

[3] Ho L L, Cavuto D J, Papavassiliou S, Zawadki A G. Adaptive and automated detection of service anomalies in transaction-oriented WANs: network analysis, algorithms, implementation, and deployment. IEEE Journal on Selected Areas in Communications, 2000, 18 (5): 744-757

[4] Li J, Liu X X, Han Z J. Research on the ARMA-based traffic prediction algorithm for wireless sensor network. Journal of Electronics and Information Technology, 2007, 29(5): 1224-1227 (in Chinese)

[5] Cadzow J A. ARMA time series modeling: an effective method. IEEE Transactions on Aerospace and Electronic Systems, 1983, AES-19(1): 49-58

[6] Versace M, Bhatt R, Hinds O, Shiffer M. Predicting the exchange traded fund DIA with a combination of genetic algorithms and neural networks. Expert Systems with Applications, 2004, 27(3): 417-425

[7] Mukherjee S, Osuna E, Girosi F. Nonlinear prediction of chaotic time series using support vector machines. In: Proceedings of the 1997 IEEE Workshop on Neural Networks for Signal Processing, 1997, 511-520

[8] Shi Z W, Han M. Support vector echo-state machine for chaotic time-series prediction. IEEE Transactions on Neural Networks, 2007, 18(2): 359-372

[9] Cao L J, Tay F E H. Support vector machine with adaptive parameters in financial time series forecasting. IEEE Transactions on Neural Networks, 2003, 14(6): 1506-1518

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg