# A Robust Data-Driven Method for Multiseasonality and Heteroscedasticity in Time Series Preprocessing

Authors: Bin Sun; Liyao Ma; Tao Shen; Renkang Geng; Yuan Zhou; Ye Tian
Journal: Wireless Communications and Mobile Computing
 (2021)

---

## Abstract

Internet of Things (IoT) is emerging, and 5G enables much more data transport from mobile and wireless sources. The data to be transmitted is too much compared to link capacity. Labelling data and transmit only useful part of the collected data or their features is a promising solution for this challenge. Abnormal data are valuable due to the need to train models and to detect anomalies when being compared to already overflowing normal data. Labelling can be done in data sources or edges to balance the load and computing between sources, edges, and centres. However, unsupervised labelling method is still a challenge preventing to implement the above solutions. Two main problems in unsupervised labelling are long-term dynamic multiseasonality and heteroscedasticity. This paper proposes a data-driven method to handle modelling and heteroscedasticity problems. The method contains the following main steps. First, raw data are preprocessed and grouped. Second, main models are built for each group. Third, models are adapted back to the original measured data to get raw residuals. Fourth, raw residuals go through deheteroscedasticity and become normalized residuals. Finally, normalized residuals are used to conduct anomaly detection. The experimental results with real-world data show that our method successfully increases receiver-operating characteristic (AUC) by about 30%.


---

## Body

## 1. Introduction

Together with rapid development of 5G, the connection requirement of wireless devices is also developing due to the eased connectivity and much shorter (in milliseconds) delay. A result is that Internet of Things (IoT) technologies are now used by more than a quarter of mainstream business compared to 13% six years ago. A great number of industry companies started to put attention on their IoT time series data, including but not limited to health care [1] and transportation [2]. While lots of mobile vehicles are connected to the IoT network as data sources [3], much more data is produced. On one aspect, it is an opportunity for machine learning-based data processing methods. On the other aspect, data transmission is now more challenging.Moving and remote data source create a challenge that it is hard to send data, especially using wireless ways, as it is still expensive to use limited wireless resource to transfer data even for 5G service providers. In some situations, if real-time moving vehicle information is needed while radio signal is limited, then wireless and wired connection may be both needed to provide support together [4, 5]. This situation is shown in Figure 1.Figure 1
Data transmission across multiple 5G wireless and wired networks among sources, edges, and centres.For this situation, one way to solve it is to label data near to sources. It not only reduces the amount of data to transfer but also balances the computing load between edges and centres [6]. One more benefit is that labelling different types of data is good for later prediction [7]. However, most solutions require labelled data to train labelling models or human expert rich experience to configurate parameters.In this work, we try to solve this problem by enhancing data preprocessing. Our previous initial feasibility experiments show promising results [8] and we complete the design here. The main contributions of this work include detailed steps of the data-driven method to handle heteroscedasticity of Internet of Things (IoT) data and comparison of possible unsupervised labelling methods as well as analysis of the reasons.The remaining content is organized as follows. First, the section introduces the problem and related definitions, together with previous research that tried to tackle this problem. Second, the proposed method is thoroughly documented in the section including steps of data preprocessing, model building, model adaptation, residual matrix construction, and anomaly detection. Third, the section describes a series of experiments using real-world data that are carried out in order to evaluate and compare the performance of the proposed method in terms of different metrics. Finally, experimental results are shown and analysed in the section, and conclusions are made in the section.


## 2. Background and Related Work

Here, we consider a system with a centre node. IoT data processing happens across the entire system [9]. It starts as early as the source application data part as shown in the updated TCP/IP architecture in Figure 2. Example source application data include camera images, video streaming, temperature, and other environmental sensed values [10]. The sensed data are then sent via possible networking routing which could be fully used for distributed processing [11], especially together with the application layer [12–14]. Physical layer choice matters as the emergency level and importance level differ among transported data which should be optimized carefully [15, 16]. When the data finally arrive at the centre, data mining algorithms could be applied [17] to analyse and conduct prediction in most cases.Figure 2
Data flow and process architecture.Regarding labelling and detection of anomalies in time series, much work has been done. Previous work can be categorized in different ways from different aspects [18]. A typical categorization includes the following categories. Probability-based methods calculate a density distribution and use some kind of thresholds to the distribution centre to label anomalies [19]. Distance-based methods set thresholds regarding how far an instance deviates from its neighbours. The measurement can be defined distances, such as in k-nearest neighbours [20], or some kind of cost of separation such as decision tree-based methods [21]. Reconstruction-based methods catch patterns and calculate the expected values of instances to get the difference, i.e., residuals, and then use residuals to conduct labelling [22, 23]. Boundary-based methods, such as support vector machine [24], provide a boundary or hyperplane to separate abnormal instances from normal ones. In addition, ensemble methods can be used to improve the accuracy and robustness of above methods [25]. For the above-mentioned methods, reconstruction-based methods give not only residuals but also comprehensive patterns and models. Thus, this work focuses on providing a preprocessing procedure to calculate and standardize residuals as the first step of reconstruction-based methods.For reconstructed residuals, as the original saved data is huge and long-term, one common problem is the variance of residuals are time-dependent, i.e., heteroscedasticity [26]. Using traffic flow as an example, the variance is high during noon time when the flow itself is high as shown in Figure 3. Vice versa, the flow and its variance are both low after midnight. This causes problems for labelling algorithms as many of them cannot distinguish high variances with anomalies.Figure 3
Overview of typical time series with heteroscedasticity.During literature review, we found two methods that try to solve the above two problems at the same time. One method is SARIMA-GARCH (Seasonal Auto-Regressive Integrated Moving Average-Generalized Auto-Regressive Conditional Heteroscedasticity) [27]. Another one is TBATS (Trigonometric Box-cox transform, ARMA errors, Trend and Seasonal component) [26]. Thus, those two methods are also tested in this work. For the final detection part, SHESD (Seasonal Hybrid Extreme Studentized Deviate test) [28] shows promising results in experiments [29–31] and is used here. It is worth mentioning that there are plenty of alternative methods while this work focuses on preprocessing.


## 3. Methodology

The proposed method includes three main steps which are preprocessing, building day-of-week (DOW) models, and solving flow-level-heteroscedasticity problem. This part describes the method in detail. The entire procedure is summarized in Figure4.Figure 4
Summary of the proposed procedure.


### 3.1. Preprocess Data

In this part, data are loaded and then divided into seven groups according to day of week.For consecutive zeros (continuous three or more zeros) which means controlled access or device malfunction, set flags and replace the instances withnull:
(1)viflag=1,ifriisinconsecutivezeros,0,otherwise.where ri is the ith measured flow rate value.Instead of using original natural daily periods, we use a new starting point. The purpose is to find a base where the starting flow rates of seasons are low and similar so that the robust fitting could work better in latter steps. It is worth mentioning that (daily) seasons may start from other time than midnight. Actually, the starting point is calculated to be around 3 am in the experiments.(2)Nseasons=NoriginNperiods,where Nseasons is the number of complete seasons, Norigin is the original number of instances (about 288×406days), and Nperiods is the number of periods (i.e., instances) per day (e.g., 288 per day for 5-minute interval data).All complete seasons are put together to construct a matrix:(3)R=s1s2⋯sis⋯sNseasons,(4)=r1,1r1,2⋯r1,is⋯r1,Nseasonsr2,1r2,2⋯r2,is⋯r2,Nseasons⋮⋮⋱⋮rip,1rip,2rip,isrip,Nseasons⋮⋮⋱⋮rNperiods,1rNperiods,2⋯rNperiods,is⋯rNperiods,Nseasons,with each season constructing a column, e.g., s1=r1,r2,⋯,rNperiodsT.Then, separate the seasons/columns into groups; here, we use the day of week of the season starting point as the criteria; thus, there are 7 groups (G1,⋯G7) with similar number of instances for each group.
(5)Gi=s7n+i∣n=0,1,2,⋯and7n+i≤Nseasons.


### 3.2. Build the Main Models

Now, seven day-of-week (DOW) models are built with the key concept of median. The building algorithm is designed in the way that it can set up several workers in parallel to improve building performance.To get a specific modelMi, a matrix R˘im is constructed by using all seasons (all columns) of Gi:
(6)R˘im=r˘1,1r˘1,2⋯r˘1,ir⋯r˘1,NMir˘2,1r˘2,2⋯r˘2,ir⋯r˘2,NMi⋮⋮⋱⋮⋮r˘ip,1r˘ip,2r˘ip,irr˘ip,NMi⋮⋮⋮⋱⋮r˘Nperiods,1r˘Nperiods,2⋯r˘Nperiods,ir⋯r˘Nperiods,NMi.NMi is the number of complete seasons for a specific model mi.Seven DOW modelsMidow (i=1,⋯,Nmodels where Nmodels is 7 in this paper) are built by applying median filters to R˘im.We can present all models as columns of a matrix:(7)M=m1m2⋯mim⋯mNmodels,(8)=r~1,1r~1,2⋯r~1,im⋯r~1,Nmodelsr~2,1r~2,2⋯r~2,im⋯r~2,Nmodels⋮⋮⋱⋮⋮r~ip,1r~ip,2r~ip,imr~ip,Nmodels⋮⋮⋮⋱⋮r~Nperiods,1r~Nperiods,2⋯r~Nperiods,im⋯r~Nperiods,Nmodels,where the im=1,2,⋯,Nmodels indicates model index and the ip=1,2,⋯,Nperiods indicates time point (period) index of day. Thus,
(9)r~ip,im=medrowipR˘im,where R˘ip,im=r˘i∈R, i.e., R˘ip,im⊂R and contains all r’s with the time point index of day ip which belongs to model Mim.


### 3.3. Adapt like Regressors

This part calculates fitted models usingM-estimation considering the above model matrix and each individual season.AnM-estimator is then computed iteratively with reweighted least squares (IRLS):
(10)βt+1=argminβ∑ip=1Nperiodswipβtεipβ2,where the scaling and addition parameters β=k,b, and residuals from the previous fit (using season is belongs to model im as an example):
(11)εβ=sis−famim,βis,im=colisR−kis,imcolimM+bis,im.Thus, the residual matrix:(12)E=ε1,1ε1,2⋯ε1,is⋯ε1,Nseasonsε2,1ε2,2⋯ε2,is⋯ε2,Nseasons⋮⋮⋱⋮εip,1εip,2εip,isεip,Nseasons⋮⋮⋱⋮εNperiods,1εNperiods,2⋯εNperiods,is⋯εNperiods,Nseasons.During the estimation, the weights are calculated as:(13)w=w1,⋯,wip,⋯,wNperiods=ψγε/cε,where c is a scaling factor:
(14)c=medabsεη,and ψ is in Huber family:
(15)ψγx=x,ifx≤γ,γsignx,ifx>γ,while η is a constant 0.675 and γ is 1.345 which correspond to regression estimator 95% efficiency. If M-estimation fails (rarely), then constrained M-estimation (CM) [32] is used (which is always working for our data). CM is proposed by Mendes and Tyler for regression and is more robust while keeping the same breakdown point (i.e., 1/2) though slower.


### 3.4. Construct the Residual Matrix

While having the adapted models, the raw residuals can be calculated directly. However, the raw residuals contain different variations on different flow levels. Thus, this part also removes flow-level-related heteroscedasticity.For adapted models, i.e.,fam,β, let us take values of adapted models and round them to integers then we get flow levels as integers l of each time point.
(16)L=fam,β,(17)=s^1s^2⋯s^is⋯s^Nseasons,(18)=l1,1l1,2⋯l1,is⋯l1,Nseasonsl2,1l2,2⋯l2,is⋯l2,Nseasons⋮⋮⋱⋮lip,1lip,2lip,islip,Nseasons⋮⋮⋱⋮lNperiods,1lNperiods,2⋯lNperiods,is⋯lNperiods,Nseasons.Be aware that the flow levels are rounded from adapted model values instead of measured. For example, suppose 9 am traffic is 85 in the DOW model, 90.3 in the adapted model, but only 10 in the measured traffic (due to an incident or so); then, the traffic flow level is 90, i.e., flow level is a adapted and generalized description which represents what the traffic should be during a similar day.Suppose the minimum and maximum integers (levels) inL are:
(19)lmin=minlip,is,lmax=maxlip,is,lip,is∈L,then we can generate a level vector
(20)L˘=lminlmin+1⋯lil⋯lmax,(21)=l˘1l˘2⋯l˘il⋯l˘Nlevels,which contains all integers from lmin to lmax and Nlevels=lmax−lmin+1 denotes the number of total flow levels.For all level items/values in adapted modelsL, adapted models do element-wise XNOR logic and we get a mask matrix A with ones indicating the time points/instances with flow levels of lil.
(22)Ail=lil⊕L¯=aia,ja∣ia=1,2,⋯,Nperiods;ja=1,2,⋯,Nseasons,(23)aia,ja=1,iflip,is=lil,null,otherwise.Let us apply this maskAil to E and take all the matched values then calculate the variance (standard deviation) for an arbirtary level null items and related calculation are ignored during this process.
(24)E˘il=A·E=εip,is∣aip,is=1,(25)vil=stdE˘il.The variances for different levels vary, thus heteroscedasticity. When putting all variances for all levels to get a variance/heteroscedasticity vector, note that residuals from neighbour levels are used when the amount of residuals is insufficient.(26)V=v1v2⋯vil⋯vNlevels.Later, all residualsE are divided by the time point’s level’s variance to get “normalized residuals.” First, for levels of each time point, i.e., lip,is, find its corresponding variance:
(27)i^lip,is=argwhereillip,is=l˘il,(28)v^ip,is=vi^lip,is.Generate a matrix of all residual’s corresponding variance:(29)V^=v^1,1v^1,2⋯v^1,is⋯v^1,Nseasonsv^2,1v^2,2⋯v^2,is⋯v^2,Nseasons⋮⋮⋱⋮v^ip,1v^ip,2v^ip,isv^ip,Nseasons⋮⋮⋱⋮v^Nperiods,1v^Nperiods,2⋯v^Nperiods,is⋯v^Nperiods,Nseasons.Normalized residuals are:(30)R′=rip,is′∣ip=1,2,⋯,Nperiods;is=1,2,⋯,Nseasons,where
(31)r′ip,is=rip,isv^ip,is.


### 3.5. Detect Using Normalized Residuals

Finally, normalized residuals are sent to detection algorithms. The entire procedure is also presented in pseudocode (Algorithm1).Algorithm 1: DOW-FLH Modelling for Data Preprocessing
1:procedure DOW-FLH (Original Time Series)2:  set flags for consecutive zeros ▹Handel Dirty Data3:for each day do4:    find the time point index (TPI) of the lowest flow5:end for6:  find TPIs’ median number as starts of daily seasons, e.g., 3 am7:for model mim in all DOW models do▹Build DOW Models8:    take all seasons related tomim to a group9:    remove flagged consecutive zeros10:    calculate median of grouped seasons as the modelmim11:end for12:for model mim in all DOW models do ▹Fit/Adapt to Get Scalings k and Additions b13:for each realted season do14:      remove flagged consecutive zeros15:      estimatek,b by robustly fitting mim to the season16:      rounding all values of the fitted model to integers as the season’s flow levels17:      get residuals as the difference between the fitted and the season18:end for19:end for20:for each flow level Standardize Residuals (FLH) do▹Standardize Residuals (FLH)21:    take all residuals for this flow level (or with neighbours if not enough)22:    calculate standard deviations (STD)23:end for24:  consider all STDs with all flow levels as the flow level heteroscedasticity (FLH)25:  divide each residual with timely corresponding STD to standardize26:for each detection algorithm do ▹Detection27:    feed the entire standardized residual time series to the algorithm28:    get algorithm-specific anomalies or anomaly scores29:end for30:return the list of anomalies or anomaly scores31:end procedure


## 3.1. Preprocess Data

In this part, data are loaded and then divided into seven groups according to day of week.For consecutive zeros (continuous three or more zeros) which means controlled access or device malfunction, set flags and replace the instances withnull:
(1)viflag=1,ifriisinconsecutivezeros,0,otherwise.where ri is the ith measured flow rate value.Instead of using original natural daily periods, we use a new starting point. The purpose is to find a base where the starting flow rates of seasons are low and similar so that the robust fitting could work better in latter steps. It is worth mentioning that (daily) seasons may start from other time than midnight. Actually, the starting point is calculated to be around 3 am in the experiments.(2)Nseasons=NoriginNperiods,where Nseasons is the number of complete seasons, Norigin is the original number of instances (about 288×406days), and Nperiods is the number of periods (i.e., instances) per day (e.g., 288 per day for 5-minute interval data).All complete seasons are put together to construct a matrix:(3)R=s1s2⋯sis⋯sNseasons,(4)=r1,1r1,2⋯r1,is⋯r1,Nseasonsr2,1r2,2⋯r2,is⋯r2,Nseasons⋮⋮⋱⋮rip,1rip,2rip,isrip,Nseasons⋮⋮⋱⋮rNperiods,1rNperiods,2⋯rNperiods,is⋯rNperiods,Nseasons,with each season constructing a column, e.g., s1=r1,r2,⋯,rNperiodsT.Then, separate the seasons/columns into groups; here, we use the day of week of the season starting point as the criteria; thus, there are 7 groups (G1,⋯G7) with similar number of instances for each group.
(5)Gi=s7n+i∣n=0,1,2,⋯and7n+i≤Nseasons.


## 3.2. Build the Main Models

Now, seven day-of-week (DOW) models are built with the key concept of median. The building algorithm is designed in the way that it can set up several workers in parallel to improve building performance.To get a specific modelMi, a matrix R˘im is constructed by using all seasons (all columns) of Gi:
(6)R˘im=r˘1,1r˘1,2⋯r˘1,ir⋯r˘1,NMir˘2,1r˘2,2⋯r˘2,ir⋯r˘2,NMi⋮⋮⋱⋮⋮r˘ip,1r˘ip,2r˘ip,irr˘ip,NMi⋮⋮⋮⋱⋮r˘Nperiods,1r˘Nperiods,2⋯r˘Nperiods,ir⋯r˘Nperiods,NMi.NMi is the number of complete seasons for a specific model mi.Seven DOW modelsMidow (i=1,⋯,Nmodels where Nmodels is 7 in this paper) are built by applying median filters to R˘im.We can present all models as columns of a matrix:(7)M=m1m2⋯mim⋯mNmodels,(8)=r~1,1r~1,2⋯r~1,im⋯r~1,Nmodelsr~2,1r~2,2⋯r~2,im⋯r~2,Nmodels⋮⋮⋱⋮⋮r~ip,1r~ip,2r~ip,imr~ip,Nmodels⋮⋮⋮⋱⋮r~Nperiods,1r~Nperiods,2⋯r~Nperiods,im⋯r~Nperiods,Nmodels,where the im=1,2,⋯,Nmodels indicates model index and the ip=1,2,⋯,Nperiods indicates time point (period) index of day. Thus,
(9)r~ip,im=medrowipR˘im,where R˘ip,im=r˘i∈R, i.e., R˘ip,im⊂R and contains all r’s with the time point index of day ip which belongs to model Mim.


## 3.3. Adapt like Regressors

This part calculates fitted models usingM-estimation considering the above model matrix and each individual season.AnM-estimator is then computed iteratively with reweighted least squares (IRLS):
(10)βt+1=argminβ∑ip=1Nperiodswipβtεipβ2,where the scaling and addition parameters β=k,b, and residuals from the previous fit (using season is belongs to model im as an example):
(11)εβ=sis−famim,βis,im=colisR−kis,imcolimM+bis,im.Thus, the residual matrix:(12)E=ε1,1ε1,2⋯ε1,is⋯ε1,Nseasonsε2,1ε2,2⋯ε2,is⋯ε2,Nseasons⋮⋮⋱⋮εip,1εip,2εip,isεip,Nseasons⋮⋮⋱⋮εNperiods,1εNperiods,2⋯εNperiods,is⋯εNperiods,Nseasons.During the estimation, the weights are calculated as:(13)w=w1,⋯,wip,⋯,wNperiods=ψγε/cε,where c is a scaling factor:
(14)c=medabsεη,and ψ is in Huber family:
(15)ψγx=x,ifx≤γ,γsignx,ifx>γ,while η is a constant 0.675 and γ is 1.345 which correspond to regression estimator 95% efficiency. If M-estimation fails (rarely), then constrained M-estimation (CM) [32] is used (which is always working for our data). CM is proposed by Mendes and Tyler for regression and is more robust while keeping the same breakdown point (i.e., 1/2) though slower.


## 3.4. Construct the Residual Matrix

While having the adapted models, the raw residuals can be calculated directly. However, the raw residuals contain different variations on different flow levels. Thus, this part also removes flow-level-related heteroscedasticity.For adapted models, i.e.,fam,β, let us take values of adapted models and round them to integers then we get flow levels as integers l of each time point.
(16)L=fam,β,(17)=s^1s^2⋯s^is⋯s^Nseasons,(18)=l1,1l1,2⋯l1,is⋯l1,Nseasonsl2,1l2,2⋯l2,is⋯l2,Nseasons⋮⋮⋱⋮lip,1lip,2lip,islip,Nseasons⋮⋮⋱⋮lNperiods,1lNperiods,2⋯lNperiods,is⋯lNperiods,Nseasons.Be aware that the flow levels are rounded from adapted model values instead of measured. For example, suppose 9 am traffic is 85 in the DOW model, 90.3 in the adapted model, but only 10 in the measured traffic (due to an incident or so); then, the traffic flow level is 90, i.e., flow level is a adapted and generalized description which represents what the traffic should be during a similar day.Suppose the minimum and maximum integers (levels) inL are:
(19)lmin=minlip,is,lmax=maxlip,is,lip,is∈L,then we can generate a level vector
(20)L˘=lminlmin+1⋯lil⋯lmax,(21)=l˘1l˘2⋯l˘il⋯l˘Nlevels,which contains all integers from lmin to lmax and Nlevels=lmax−lmin+1 denotes the number of total flow levels.For all level items/values in adapted modelsL, adapted models do element-wise XNOR logic and we get a mask matrix A with ones indicating the time points/instances with flow levels of lil.
(22)Ail=lil⊕L¯=aia,ja∣ia=1,2,⋯,Nperiods;ja=1,2,⋯,Nseasons,(23)aia,ja=1,iflip,is=lil,null,otherwise.Let us apply this maskAil to E and take all the matched values then calculate the variance (standard deviation) for an arbirtary level null items and related calculation are ignored during this process.
(24)E˘il=A·E=εip,is∣aip,is=1,(25)vil=stdE˘il.The variances for different levels vary, thus heteroscedasticity. When putting all variances for all levels to get a variance/heteroscedasticity vector, note that residuals from neighbour levels are used when the amount of residuals is insufficient.(26)V=v1v2⋯vil⋯vNlevels.Later, all residualsE are divided by the time point’s level’s variance to get “normalized residuals.” First, for levels of each time point, i.e., lip,is, find its corresponding variance:
(27)i^lip,is=argwhereillip,is=l˘il,(28)v^ip,is=vi^lip,is.Generate a matrix of all residual’s corresponding variance:(29)V^=v^1,1v^1,2⋯v^1,is⋯v^1,Nseasonsv^2,1v^2,2⋯v^2,is⋯v^2,Nseasons⋮⋮⋱⋮v^ip,1v^ip,2v^ip,isv^ip,Nseasons⋮⋮⋱⋮v^Nperiods,1v^Nperiods,2⋯v^Nperiods,is⋯v^Nperiods,Nseasons.Normalized residuals are:(30)R′=rip,is′∣ip=1,2,⋯,Nperiods;is=1,2,⋯,Nseasons,where
(31)r′ip,is=rip,isv^ip,is.


## 3.5. Detect Using Normalized Residuals

Finally, normalized residuals are sent to detection algorithms. The entire procedure is also presented in pseudocode (Algorithm1).Algorithm 1: DOW-FLH Modelling for Data Preprocessing
1:procedure DOW-FLH (Original Time Series)2:  set flags for consecutive zeros ▹Handel Dirty Data3:for each day do4:    find the time point index (TPI) of the lowest flow5:end for6:  find TPIs’ median number as starts of daily seasons, e.g., 3 am7:for model mim in all DOW models do▹Build DOW Models8:    take all seasons related tomim to a group9:    remove flagged consecutive zeros10:    calculate median of grouped seasons as the modelmim11:end for12:for model mim in all DOW models do ▹Fit/Adapt to Get Scalings k and Additions b13:for each realted season do14:      remove flagged consecutive zeros15:      estimatek,b by robustly fitting mim to the season16:      rounding all values of the fitted model to integers as the season’s flow levels17:      get residuals as the difference between the fitted and the season18:end for19:end for20:for each flow level Standardize Residuals (FLH) do▹Standardize Residuals (FLH)21:    take all residuals for this flow level (or with neighbours if not enough)22:    calculate standard deviations (STD)23:end for24:  consider all STDs with all flow levels as the flow level heteroscedasticity (FLH)25:  divide each residual with timely corresponding STD to standardize26:for each detection algorithm do ▹Detection27:    feed the entire standardized residual time series to the algorithm28:    get algorithm-specific anomalies or anomaly scores29:end for30:return the list of anomalies or anomaly scores31:end procedure


## 4. Experiments

This section describes data, practical procedure, and the way we conduct experiments.


### 4.1. Data Specification

The one-year long real-world data are collected from a highway. Ground truth anomaly (incidents) labels are generated by using the extended system mentioned in [33]. The data are imputed using the method from [34] before any processing. One device sends a monitored flow record at five-minute intervals. Each record contains some traffic statistics such as flow rate and average speed. This road carries undersaturated flow except in holidays’ noons, where is 15 min average?


### 4.2. Experimental Setup

The experiments are done in a desktop computer with AMD Ryzen 5-3600 (6 Cores, 3600 MHz) and 16 GB DDR4 memory. To be fair, we only implement our method; other algorithms are taken from public domain such as GitHub.Our implementation is done in the R programming environment version 3.4.3 with RStudio 1.3.1056, AnomalyDetection 1.0, forecast 8.2, feather 0.3.3 as well as the Python programming environment version 3.6.7/3.6.9 with library arch 4.8.1, statsmodels 0.9.0, feather-format 0.4.0/0.4.1, numpy 1.16.0/1.19.4, pandas 0.23.4/1.1.4, scikit-learn 0.19.2/0.23.2, scipy 1.2.2/1.5.4, ipykernel 5.3.4, and ipython 7.16.1.SHESD was originally implemented to give only binary results so we modified it by addingtestresult−criticalvalue to get anomaly/outlier scores. Also, as the max allowed anomaly (outliers) ratio is 50%, we mark all nontested ones the same score as the lowest score.


### 4.3. Evaluation Measurement and Metrics

Receiver-operating characteristic (ROC) is used as the main evaluation metric as it provides an accurate and visualized way to present detecting results. One important value from ROC is area under curve (AUC) which is also known as A′ (“a-prime”), or concordance-statistic (c-statistic). It is a measure of goodness of fit that is often used for binary classification modelling results evaluation; therefore, we use it here.


## 4.1. Data Specification

The one-year long real-world data are collected from a highway. Ground truth anomaly (incidents) labels are generated by using the extended system mentioned in [33]. The data are imputed using the method from [34] before any processing. One device sends a monitored flow record at five-minute intervals. Each record contains some traffic statistics such as flow rate and average speed. This road carries undersaturated flow except in holidays’ noons, where is 15 min average?


## 4.2. Experimental Setup

The experiments are done in a desktop computer with AMD Ryzen 5-3600 (6 Cores, 3600 MHz) and 16 GB DDR4 memory. To be fair, we only implement our method; other algorithms are taken from public domain such as GitHub.Our implementation is done in the R programming environment version 3.4.3 with RStudio 1.3.1056, AnomalyDetection 1.0, forecast 8.2, feather 0.3.3 as well as the Python programming environment version 3.6.7/3.6.9 with library arch 4.8.1, statsmodels 0.9.0, feather-format 0.4.0/0.4.1, numpy 1.16.0/1.19.4, pandas 0.23.4/1.1.4, scikit-learn 0.19.2/0.23.2, scipy 1.2.2/1.5.4, ipykernel 5.3.4, and ipython 7.16.1.SHESD was originally implemented to give only binary results so we modified it by addingtestresult−criticalvalue to get anomaly/outlier scores. Also, as the max allowed anomaly (outliers) ratio is 50%, we mark all nontested ones the same score as the lowest score.


## 4.3. Evaluation Measurement and Metrics

Receiver-operating characteristic (ROC) is used as the main evaluation metric as it provides an accurate and visualized way to present detecting results. One important value from ROC is area under curve (AUC) which is also known as A′ (“a-prime”), or concordance-statistic (c-statistic). It is a measure of goodness of fit that is often used for binary classification modelling results evaluation; therefore, we use it here.


## 5. Results and Analysis

As shown in Figure5, our DOW and DOW-FLH methods are superior with regard to AUC. DOW with and without FLH performs similar considering AUC of 0.693 from both algorithms which are 26.9% better coverage than other algorithms on average (AUC 0.546). What is more, DOW-FLH is preferred for less false positives on the optimal cut-off point compared to DOW without FLH due to the data sensitivity to false positive. May move below to analysis? For unbalanced datasets such as traffic flows, this behaviour gives positive influence. The reason is that some false-negative instances introduce only minor issues for true-negative ones as negative instances are majority while the same amount false-positive instances impact true anomalies (incidents) much more.Figure 5
The proposed DOW and DOW-FLH gives similar AUC and around 27% bigger coverage than the others on average. Besides, DOW-FLH produces only half false positives compared to DOW without FLH.We analysed the detection ratio and AUCs for different situations and found some interesting results. For device malfunction incidents, most algorithms cannot notice it as shown in Figure6. The possible reason is that other algorithms are tracking no-flow situation without considering normal situation. Note that good seasonal modelling (DOW) should work with suitable variance handling methods, as inappropriate variation handling (i.e., GARCH) may otherwise reduce the effectiveness.Figure 6
Only DOW-based algorithms can detect device malfunction efficiently while others are close to no discrimination line.Figure7 shows level to residual characteristics diagnostics. The mean of residuals (blue line) is mostly under 2 but increases rapidly to be about 5 when the flow level is greater than 150. This is due to the fact that extreme levels (greater than 150) occur only during few big holidays, so this scenario is hard to be caught by models. The standard deviations (green line) is mainly increasing which represents one key problem, i.e., heteroscedasticity. The purple line represents the number of instances per level, and it becomes very small for extreme scenarios in both directions of x-axis. The number of span is used to include neighbour levels when one level’s corresponding instances are too few to calculate reasonable statistics. In summary, it can be seen that the relation mapping from levels to residual characteristics are nonlinear. This explains why the proposed data-driven algorithms perform better.Figure 7
Flow level with corresponding instance numbers and mean, STD, and span of residuals (instance number per level is scaled down to be shown in this figure).DOW successfully modelled patterns and FLH successfully suppressed heteroscedasticity for normal data compared to others as the residuals are shown in Figure8. Other algorithms, when being compared to DOW-FLH, cannot distinguish data with vs. without abnormalities, such as shown in Figure 9. This could be an advantage for GARCH-based methods when tracing rapid change in (nonseasonal) time series with heteroscedasticity, but it becomes an disadvantage and hides possible abnormal data instances here. The problem with TBATS and SARIMA is that they could not successfully model the patterns and produces residuals with much noise which leads to low signal-noise ratio as shown in Figure 10.Figure 8
The left one shows DOW-FLH residuals of time series data without and with abnormal data and the right one shows overall residual density for the entire dataset. DOW-FLH successfully suppressed heteroscedasticity for normal data.Figure 9
DOW-GARCH suppresses not only extreme values and heteroscedasticity for data with all normal instances but also for data with abnormal instances, when being compared to DOW-FLH (Figure9) (same data range).Figure 10
TBATS-based algorithms produce residuals with much noise when being compared to DOW-FLH (Figure9) (same x-axis time range). (S)ARIMA-based algorithms give similar results.Previous work has shown that ARIMA and GARCH cannot be adapted to seasonality with many periods such as here 288 periods per season. Instead, they will adapt to local trend or rapid change add plots; therefore, they are not suitable to detect anomalies lasting beyond their detection abilities. This characteristic could be an advantage when quick predicting traffic for short-term time is needed.


## 6. Conclusion and Future Work

The experiment results show that the proposed DOW algorithm is good at matching multiseasonality time series patterns, and FLH can solve heteroscedasticity problem. DOW-FLH-modelled residuals can be used for labelling anomalies; then, the chosen data can be sent to either edges or centres for further process.As discussed above, the proposed DOW-FLH in this work is good at modelling and labelling multiseasonal IoT time series for the edge-centre structure. However, other compared algorithms, including SARIMA- and TBATS-based ones, are more mature and may be good at local trend prediction. Also, edge computing can engage crowdsourcing and related active learning [35] to make full use of advantages provided by edge-centre structure.This point can be further tested in later research.Labelling can be treated as a classification question, and many new algorithms can work on this task. Especially, recent development regarding classification using belief theory is showing promising results [36], and it is good for multisource scenarios in edge-centre computing. Thus, this might be a good enhancement for our current work, and we look forward to investigate more about it in the future work.In summary, the proposed DOW-FLH method performs well during experiments using multiseasonal IoT time series and should be considered to use when labelling is needed in edge-centre computing structure.



---

*Source: 6692390-2021-08-16.xml*