Improving Time Series Forecasting by Discovering Frequent Episodes in Sequences

Improving Time Series Forecasting by Discovering Frequent Episodes in Sequences Francisco Mart´ınez-´ Alvarez1, Alicia Troncoso1,andJos´e C. Riquelme2 1Area of Computer Science, Pablo de Olavide University of Seville, Spain {fmaralv,ali}@upo.es 2Department of Computer Science, University of Seville, Spain [email protected] Abstract. This work aims to improve an existing time series forecasting algorithm –LBF– by the application of frequent episodes techniques as a complementary step to the model. When real-world time series are forecasted, there exist many samples whose values may be specially unexpected. By the combination of frequent episodes and the LBF algorithm, the new procedure does not make better predictions over these outliers but, on the contrary, it is able to predict the apparition of such atypical samples with a great accuracy. In short, this work shows how to detect the occurrence of anomalous samples in time series improving, thus, the general forecasting scheme. Moreover, this hybrid approach has been successfully tested on electricity-related time series. Keywords: Time series, forecasting, outliers. 1 Introduction This work provides a new methodology to forecast time series and, in addition, to predict the apparition of outliers. The analysis of temporal data and the forecast of future values of time series are among the most important problems that data analysts face in many fields, ranging from finance and economics, to production operations management or telecommunications. A forecast is a prediction of some future events. The proposed approach is specifically framed in electricity prices time series forecasting, which is a difficult task due to the nonconstant mean and variance and significant outliers typically present in these series. Thus, the combination of two different techniques are proposed to fulfill this goal. The first one is a general-purpose forecasting algorithm introduced in [9], called LBF. The authors obtained a previous labeling of the elements forming the time series by means of clustering techniques. The forecasting process was performed by using just the information provided by the clustering. Thus, the values of the elements in datasets were discretized and, as a result, the sequence of real values was transformed in a sequence of discrete values or labels. These labels were used to predict the future behavior of the time series, avoiding the N. Adams et al. (Eds.): IDA 2009, LNCS 5772, pp. 357–368, 2009. c Springer-Verlag Berlin Heidelberg 2009 358 F. Mart´ınez-´ Alvarez, A. Troncoso, and J.C. Riquelme use of the real values until the last step of the process. The results returned by the algorithm, however, were not labels but the real values. The algorithm introduced in [10] is inserted in the general scheme of forecasting with the aim of dealing with the presence of outliers. Concretely, they proposed an algorithm called Q-epiMiner that was, in fact, an improvement of the well-known serial episodes [8]. The main achievement of the Q-epiMiner algorithm was to characterize sequences of similar behavior over all occurrences, as well as providing a tree structure to organize these sequences. Therefore, the discovery of frequent episodes is used in order to determine possible candidates to be outliers when using the LBF algorithm to forecast time series. The general outline of the new proposed approach is shown in Fig. 1. Discretization (discrete values) Time series (real values) Frequent episodes (special sequences) Forecasting process Technique in [9] Technique in [9]Technique in [10] Fig. 1. General outline of the proposed methodology The rest of the paper is organized as follows. The latest works related to time series forecasting and discovery of outliers are gathered in Section 2. Section 3 presents a brief explanation of the two existing and used algorithms in which the new approach is based on. Section 4 introduces the proposed methodology, showing how the two existing techniques are combined in order to improve the forecasting process. Section 5 shows the results obtained for the electric energy market of Spain for the year 2006, including measures of the quality of them. Finally, Section 6 summarizes the main conclusions achieved and provides clues for future work. 2 Related Work Regarding to time series forecasting, two mixed models were proposed to obtain the forecasting of the prices for two different prediction horizons in [3]. The first one forecasted electricity prices for each of the 24 hours of the next day using ARIMA models, while the second model computed the predictions by using Bayesian information criteria. A modification of the nearest neighbors methodology was proposed in [13]. To be precise, the approach presented a simple technique to forecast next-day electricity market prices based on the weighted nearest neighbors methodology. Li et al. proposed a forecasting system immersed in a grid environment in [6]. In this paper, a fuzzy inference system –adopted due to its transparency and interpretability– and time series methods were proposed for day-ahead electricity price forecasting. The authors in [12] proposed an artificial neural network-based approach to forecast the energy price in the Spanish market. Thus, a novel training method Improving Time Series Forecasting by Discovering Frequent Episodes 359 was presented and applied to the multilayer perceptron in order to improve the forecasting process. The apparition of outliers in time series has also been widely discussed in the literature. Thus, the authors in [1] proposed a technique to detect outliers in data whose generation was difficult to model. They assumed that the correlation among data close in time is higher than those farther apart. Another method for detecting outliers was proposed in [7], in which the authors considered two different sources of outliers –additive and innovation– in autorregressive moving-average time series. Concretely, they proposed the application of two different procedures associated to each source simultaneous and coherently. The occurrence of spike prices (prices significantly higher than the expected value, i. e., outliers) is an usual feature associated with electricity prices-related time series. With the aim of dealing with this peculiarity, the authors in [14] proposed a data mining framework based on both support vector machines and probability classifiers. The work described in [5] searched for patterns in electricity prices data in order to verify how the outliers may modify the behavior of such prices. To fulfill this goal, they used Box and Jenkins models, Discrete Fourier Transform series smoothing and GARCH approaches. Also, the work in [4] discussed the use of the fractal theory to forecast the electricity price time series. For this purpose, a forecasting model based on improved fractal was built and solved to forecast short–term electricity price time series. 3 Fundamentals The proposed approach is a combination of two existing techniques. Thus, this section provides the mathematical fundamentals underlying to both LBF (Subsection 3.1) and Q-epiMiner algorithms (Subsection 3.2). A more detailed explanation can be found in [9] and [10], respectively. 3.1 Time Series Forecasting: The LBF Algorithm The LBF algorithm was initially presented in [9]. Given the hourly prices recorded in the past, up to day d, the forecasting problem aims to predict the 24 hourly prices corresponding to day d+1. Let Pi∈R24 be a vector composed of the 24 hourly energy prices corresponding to a certain day i Pi=[p1,p 2,...,p 24].(1) Let Libe the label of the prices of the day iobtained as a previous step to the forecasting by using a clustering technique. Let Si Wthe subsequence of labels of the prices of the Wconsecutive days, from day ibackward, as follows: Si W=[Li−W+1,L i−W+2,...,L i−1,L i](2) where the length of the window, W, is a parameter to be determined. 360 F. Mart´ınez-´ Alvarez, A. Troncoso, and J.C. Riquelme The LBF algorithm first searches for the subsequences of labels which are exactly equals to Sd Win the data base, providing the equal subsequences set, ES, defined by the equation, ES =set of indexes jsuch that Sj W=Sd W(3) In case of not finding any subsequence in data base equal to Sd W, the procedure searches for the subsequences of labels which are exactly equals to Sd W−1.That is, the length of the window composed of the subsequence of labels is decreased. According to the LBF approach, the 24 hourly prices of day d+1 are predicted by averaged the prices of the days succeeding those in ES.Thatis, Pd+1 =1 size(ES)· j∈ES Pj+1 (4) where size(ES) is the number of elements belonging to the set ES.Afterwards, LBF algorithm outputs need to be de-normalized to generate the desired forecasted values. When the horizon of prediction is greater than one day, the following tasks have to be carried out. First of all, the real values of the predicted sample are linked to the whole dataset. Second, the clustering process is repeated with the enlarged dataset and, finally, the window size is re-calculated and the prediction step is performed. 3.2 Frequent Episodes in Sequences: The Q-epiMiner The algorithm introduced in [10] analyzes data events or sequences in order to find episodes. Formally, a sequence is an ordered list of events, where an event is identified by the pair ev =< date, eventT ype >. The ordered occurrence of events is called serial episode and represented by: E=[ev1,ev 2, ..., evn]=[<d 1,t 1>, < d2,t 2>, ..., < dn,t n>](5) where nis the number of events forming an episode. Thus, the algorithm is able to handle with three time constraints, provided as input data: the minimum time span between two events or gapmin,themaximum time span between two events or gapmax and the maximum time span between the beginning and the end of an episode or windowSize. In this way, an episode has to simultaneously satisfy three constraints: the gapmin has to be greater or equal to a given threshold, the gapmax hastobelesserorequaltoagiven threshold and the windowSize has to be lesser or equal to a given threshold. The thresholds are set depending on the requirements of the application. The rules used in the algorithm are computed where the antecedent is a serial episode and the consequent contains only one event type. Then, a list of the Improving Time Series Forecasting by Discovering Frequent Episodes 361 positions in the data sequence LEof a particular episode is built. Concretely, the list contains the time stamps associated with the events comprised in the episode and they are sorted by increasing values. The next step consists in evaluating the whole sequence of events and the LE. From this evaluation, the algorithm provides a set of tuples <ev,L ev >,where Lev is the list of locations of the event occurrences and ev is an event of the episode E. Finally, the standard prefix-based strategy [2] is used for the overall enumeration since it fits well with both episodes extraction and the use of the sorted lists. In other words, an episode is used as a prefix and expanded in order to obtain new episodes. 4 Methodology This section explains the methodology proposed to improve the forecasting process provided by the LBF algorithm. The discovery of frequent episodes is included in the aforementioned algorithm as a crucial step for outliers detection. Thus, the proposed methodology is divided into two phases clearly differentiated. First of all, the LBF algorithm is trained with the datasets under study. Second, the predictions with the highest error rates made during the training are analyzed by means of frequent episodes techniques. From this analysis, the days likely to be outliers will be determined and not considered in the prediction process. 4.1 Combining the LBF Algorithm with Frequent Episodes The value of two parameters have to be determined in the LBF process, K and W. With regard to the number of clusters (K), the proposed approach acts exactly equal to what was proposed in [9] (see Section II.C). However, it is important to remark that the length of the window (W) is slightly different calculated from how it was proposed in the original paper. In practice, Wis calculated by means of cross-validation. Concretely, the n−fold cross-validation is used in this work to obtain the optimal value of W.Inn−fold cross–validation, the original sample is partitioned into nsubsets. Of the nsubsets, a single subset is retained as the validation data for testing the model, and the remaining n−1 subsets are used as training data. The cross-validation process is then repeated ntimes (the folds), with each of the nsubsets used exactly once as the validation data. The nresults from the folds are then averaged (even if some authors prefer a combination of them) to produce a single estimation. The advantage of this method over repeated random sub-sampling is that all observations are used for both training and validation, and each observation is used for validation exactly once. Twelve folds have been created in this work (n= 12) for all datasets, where each fold represents a month. Consequently, the training set consists of one year. 362 F. Mart´ınez-´ Alvarez, A. Troncoso, and J.C. Riquelme The 12−fold cross–validation is then evaluated. The forecasting errors are calculated in every fold by varying the length of W. These monthly errors are denoted by emonth{W=j}for j=1...W max,whereWmax = 10 –no longer sequences were found in datasets. Then, the average errors are calculated for each window size as follows, ei=n i=1 emonth{W=i} n(6) where n=12andmonth ={Jan,...,Dec}. The Wselected is the one that minimizes the average error corresponding to the 24 folds (months) evaluated. W=argmin{ei}with i=1, ..., Wmax (7) This modification has a simple justification. With the former methodology (see Section II.D in [9]), only one test set was evaluated and, consequently, the number of episodes found may be limited and not conclusive. However, with the application of n−fold cross-validation, the number of sequences is increased n times providing thus ntraining sets instead of one. It is now –just after the training step and before the prediction process– that the discovery of frequent episodes plays an important role. Hence, the apparition of anomalous days is intended to be predicted. Once the n−fold cross–validation is applied, the number of sequences generated can be calculated as: #S(W, FL)=nFL−W+1 (8) where nis the number of folds, FL the average fold length and Wthe length of the sequence (or window) considered. Note that the maximum number of dissimilar sequences that can be generated is bounded by: Nmax(K, W)=KW(9) which will be a number typically much higher than the sequences found in the training sets (Nmax >> #S(W, FL)). The episodes aimed to be found are those which generate a prediction error greater than the average error in the cross–validation process. For this reason, a set of events that satisfies ej>min{ei},fori=1, ..., Wmax and j=1, ..., nFL, is constructed. This set, CS, gathers all the candidates events to be preceded by an episode precursor of outliers. Nevertheless, not all these candidates have the same probability to be outliers since the associated errors range from values near to the mean error (these candidates should be finally discarded by the approach) to values significantly high. For this reason, each candidate is co-labeled by using clustering techniques, concretely, the K-means algorithm. The decision on how many clusters have to be created is always an open question and many indices might be used. However, Improving Time Series Forecasting by Discovering Frequent Episodes 363 Data acquisition Normalization Clustering PREDICTION Select K (Silhouette) Denormalization More days? End Yes No Labeleddata Select W (12-fold CV) DISCOVERY OF EPISODES Episodes Predicted day insertion Fig. 2. Illustration of the proposed methodology Construct CS Clustering CS Find episodes Candidates Prediction DISCOVERY OF EPISODES {Cl,C m,C u} Clustering Set W Search for W In historical data Calculate sample-ahead Denormalization PREDICTION Discovery of episodes W matches episodes ? OUTLIER YES NO Fig. 3. Detail of both discovery of episodes and prediction steps it is worthless to have a large number of clusters and for this reason only three conceptual classes will be created: lower (Cl), central (Cc) and upper (Cu) errors. An a priori reasoning reveals that those candidates belonging to Cumust be more probable to be preceded by sequences generators of outliers than the candidates in Clor Cc. Results corresponding to each cluster of data will be separately analyzed in Section 5. The next step consists in computing the episodes (concrete sequences of labels) occurred before the candidates in order to determine the apparition of an outlier. Hence, the approach has to decide which sequences preceding the candidates are the episodes causing errors greater than the expected average error. Then, the sequences that only appear before the candidates are considered frequent episodes and therefore preceding outliers. Fig. 2 illustrates the whole process of prediction when the frequent episodes are included in the LBF algorithm. In addition, the steps corresponding to discovery of episodes and prediction are further detailed in Fig. 3. 364 F. Mart´ınez-´ Alvarez, A. Troncoso, and J.C. Riquelme 4.2 Parameters of Quality The parameters used in order to measure the accuracy of the approach are now introduced. Note that in subsequent equations, true positives or TP is the number of candidates that indeed were preceded by episodes that caused errors greater than the average; true negatives or TN is the number of sequences found before a candidate that was properly discarded; false positives or FP is the number of candidates whose preceding episodes were erroneously considered to be causing of errors greater than the average and, finally, false negatives or FN is the number days not considered candidates and eventually preceded by episodes causing errors greater than the average. According to these definitions, the sensitivity is the probability to detect a frequent episode as precursor of outliers. Its formula is defined as follows: Sensitivity =TP TP +FN (10) Other parameter is the specificity which is the ratio of outliers candidates properly discarded by the approach. The mathematical expression is: Specificity =TN TN +FP (11) The positive predictive value (PPV) is the probability that a detected outlier is indeed a real one. Its formula is: PPV =TP TP +FP (12) Finally, the negative predictive value (NPV) is the probability that a discarded candidate to be outlier was not indeed a real one. Its formula is: NPV =TN TN +FN (13) 5Results In order to prove that the LBF works properly over different datasets, the authors in [9] considered several public electricity prices time series. The new approach is applied on the Spanish electricity price time series (OMEL) [11]. This section is structured as follows. First, the training of the LBF is presented, obtaining thus the adequate values for the parameters Kand W.The results provided in this step are, then, analyzed by means of frequent episodes techniques intending to find those patterns in the historical data that perform the worst predictions. From this analysis, some days will become candidates to have an anomalous behavior and, consequently, have a higher error Improving Time Series Forecasting by Discovering Frequent Episodes 365 prediction. Finally, the validity of considering these days candidates to be outliers is discussed. 5.1 Discovering Frequent Episodes in Time Series The forecasting process is applied for the year 2006, with a historical data of a length of one year and with a horizon of prediction of one month. Given this situation, every time a month is forecasted the training set changes. When January 2006 is forecasted, the training set comprises the whole year 2005. However, when February 2006 is forecasted the historical data ranges from February 2005 and January 2006, and so on. The results of January 2006 are now described since the explanation of the remaining eleven months is analogous. As for the LBF, the number of clusters to be generated as well as the length of the window to be searched for are calculated according to the methodology presented in Section 4.1. Thus, this pair of parameters are equal to: (K, W)←(4,5). Consequently, the number of sequences generated during the training step is #S(W, FL) = 361. However, many sequences appeared repeatedly and the final number of different sequences were 43. As the maximum number of possible sequences is Nmax(K, W) = 1024, the aforementioned number of sequences represent the 4.19% of the potential. With regard to the Q-epiMiner, the parameters are set to gapmax =1, gapmin =1andwindowSize =Win order to adapt its application to the particular problem tackled in this work. Also note, that the events are the labels generated during the LBF process, the date is the day associated with such label and the type of event is the curve of prices associated to this day. The CS can be now constructed. For this purpose, the eifrom (5) have to be considered since the candidates are those days belonging to the training set that obtained an error greater than ei. The value of the mean error, calculated according to the methodology in Section 4.1 is e=2.23%. From the 365 days comprising the training set, 131 had an error greater than 2.23% so the constructed CS contains 131 candidates. The candidates have to be classified by means of K-means, with K=3as discussed in Section 4.1. The obtained cutoff values were 3.78% and 5.69% for dividing classes Cl–Cmand Cm–Cu, respectively. From these cutoffs, the candidates were classified as follows: 98 ∈Cl,25∈Cmand 8 ∈Cu. Once the candidates are selected and classified, the sequences that generated them are evaluated. From the candidates in Cl, 7 different sequences were found (#Sl= 7); from the candidates in Cm,4(#Sm= 4) and from the candidates in Cu,2(#Su= 2). This fact involves that from the 43 sequences found in the training set, only 13 caused errors greater than the average. Finally, the number of episodes causing outliers are determined. From the sequences Cl, just one appeared exclusively as a precursor of an outliers. With reference to sequences in Cm, three out of four. And both two sequences in Cu were exclusive. Table 1 summarizes the results for the twelve months of the year 2006.