Discovery of motifs to forecast outlier occurrence in time series

Discovery of motifs to forecast outlier occurrence in time series F. Martínez–Álvarez , A. Troncoso, J.C. Riquelme , J.S. Aguilar–Ruiz Keywords: Time series forecasting Pattern recognition Motifs Outliers abstract The forecasting process of real-world time series has to deal with especially unexpected values, commonly known as outliers. Outliers in time series can lead to unreliable modeling and poor forecasts. Therefore, the identification of future outlier occurrence is an essential task in time series analysis to reduce the average forecasting error. The main goal of this work is to predict the occurrence of outliers in time series, based on the discovery of motifs. In this sense, motifs will be those pattern sequences preceding certain data marked as anomalous by the proposed metaheuristic in a training set. Once the motifs are discovered, if data to be predicted are preceded by any of them, such data are identified as outliers, and treated separately from the rest of regular data. The forecasting of outlier occurrence has been added as an additional step in an existing time series forecasting algorithm (PSF), which was based on pattern sequence similarities. Robust statistical methods have been used to evaluate the accuracy of the proposed approach regarding the forecasting of both occurrence of outliers and their corresponding values. Finally, the methodology has been tested on six electricity-related time series, in which most of the outliers were properly found and forecasted. 1. Introduction This work proposes a new strategy to predict the occurrence of outlying data in time series, as well as providing accurate forecasts for them. It is worth highlighting that the goal of this methodology is to forecast their appearance, instead of detecting them in an already known set of values, which is a common goal in robust statistics (Maronna et al., 2007). The majority of robust statistical techniques perform a posteriori detection, that is, they determine whether a datum is an outlier or not, but once it has already occurred. However, a comparison with these techniques will be of the utmost importance in order to evaluate the accuracy of the proposed metaheuristic. A general-purpose forecasting algorithm, called PSF, was presented in Martínez-Álvarez et al. (in press). Its main feature lied in performing a discretization of the time series by means of certain clustering technique. Then, it only used the generated labels to make predictions. Based on that previous discretization, this work attempts to discover pattern sequences (henceforth called motifs) in the historical data to forecast the occurrence of outliers and their associated values. This novel methodology is inserted in the general scheme of PSF. The prediction of outliers plays an important role as wrong models and poor forecasts are obtained when ignoring outliers. This work presents a metaheuristic to discover motifs in time series and,then,ifdatatobepredictedareprecededbyanyofthesediscovered motifs, consider these data as outlier. The motifs are determined during the training phase as those pattern sequences that precede data with remarkable forecasting error. Thus, the existing PSF algorithm is modified by adding a new motif extraction step. The enhanced version is capable of predicting the appearance of such outliers with great reliability when the motifs extraction step is added. In fact, the approach has been successfully tested on six real-world electricity-related time series, in particular, on energy prices and demand of three different markets, reaching sensitivity values greater than 82%, and specificity values greater than 95%. Furthermore, results about the effect of outliers on the average forecasting errors are reported for all the six time series, exhibiting remarkable forecasting error reduction. Despite the vast variety of works related to outliers detection and motifs discovery in time series, there is no approach in time series in order to forecast the occurrence of outliers, to the authors’ knowledge. The remaining of the paper is organized as follows. A review of the most recently published works regarding energy time series forecasting, motifs discovery and outliers detection can be found in Section 2. Section 3provides formal description for sensitive terms, and presents a brief explanation of the original algorithm. As for Section 4, it introduces the proposed methodology, showing how to insert the outlier occurrence forecasting in the original algorithm’s general scheme. The results obtained for the six electricity prices and demand time series are reported and discussedin Section 5. Finally, Section 6summarizes the main conclusions achieved. 2. Related work This section provides useful and recent references about the three main topics involved in this paper: Energy time series forecasting, motifs discovery in temporal data and robust statistical methods to detect outliers. For the sake of clarity, these topics have been separated in three different sections. 2.1. Energy time series forecasting The interest of analyzing electricity price time series resides in the progressive deregulation of electric power markets. Furthermore, electricity price time series possess certain features that turn the prediction into a difficult task: non-constant mean/variance and frequently outlier occurrences. For this reason, electricity-pro- ducer companies want optimized bidding strategies as well as needing assessment about the risk of trusting forecasts (Plazas et al., 2005). On the other hand, the process of forecasting the quantity of electricity required for a specific geographical area during a time period is called load forecasting or demand forecasting. This process is key since current technology allows to store a limited amount of electricity in batteries. Therefore, the demand forecasting plays an important role for electricity power suppliers because both excess and insufficient energy production may lead to increased costs and a significant reduction of profits. The pursuit of accurate forecasting in electricity price time series has motivated research works by many authors (Aggarwal et al., 2009). Thus, the use of mixed models was proposed in García- Martos et al. (2007) to forecast prices for different horizons of prediction. Also remarkable was the work introduced in Troncoso et al. (2007) that, by means of weighted nearest neighbors methodology, forecasted next-day electricity prices. Also, the use of an artificial neural network to fulfil the same goal can be found in Pino et al. (2008). Even the use of classical autoregressive models has been recently used to forecast prices in several markets (Weron and Misiorek, 2008). On the contrary, the authors in Troncoso et al. (2004) proposed a weighted nearest neighbors-based methodology to forecast electricity demand. The proposed approach was tested over the nextday Spanish load forecasting. Also, the authors in El-Telbany and El-Karmi (2008) forecasted the Jordanian electricity demand with an artificial neural network, which was trained by means of particle swarm optimization techniques. This market was also studied in Badran et al. (2008) but, this time, the authors preferred to concentrate on short and medium-term load forecasting by using regression models. Finally, Wang and Wang (2008) proposed a new prediction approach based on support vector machines (SVM) techniques with a previous selection of features from data sets by using an evolutionary method. The discovery of outliers in electricity price time series has also been widely discussed in literature. Thus, the authors in Lu et al. (2005) proposed a model based on the analysis of several variables (among which the electricity demand is highlighted) by means of Bayesian classification (BC) and similarity searching techniques. A hybrid methodology that combined SVMs and BC was developed in Wu et al. (2006) to classify both spikes and normal electricity prices. Alternatively, a data mining framework based on SVM and probability classifiers was described in Zhao et al. (2007) with the aim of forecasting spikes in prices accurately. By contrast, the prediction of peaks in electricity demand was addressed in Saini (2008). This work forecasted demand peaks up to seven days ahead using feed forward neural network and adaptive backpropagation learning methods. In the work introduced in Ismail et al. (2009), the authors developed a rule-based method that combined regression models and fuzzy systems to analyze daily electricity peak load demands in Malaysia. Also, the authors in Hyndman and Fan (2010) described a semi-parametric additive model to discover relationships between the demand and exogen variables. The approach was applied to long-term peaks for the South Australian market. 2.2. Motifs discovery The discovery of motifs in continuous data, also known as functional data in many works (Valderrama, 2008), was originally formalized in Lin et al. (2002), in which the authors introduced several algorithms to mine motifs in time series, among which the k-motif algorithm highlights. However, the main drawback of this algorithm is its dependence on a pre-fixed pattern length. Later, the authors in Tang and Liao (2008) proposed a modified version that improved, precisely, this feature. Furthermore, they generated original patterns by considering the discovered motifs. The detection of on-line motifs in continuous data has also been addressed. Particularly, a new methodology to detect on-line motifs in time series by combining probabilistic models and polynomial least-squares approximations was proposed in Fuchs et al. (2009). This topic was also studied in Mueen and Keogh (2010) in which the authors found and maintained time series motifs from robotics, online compression and wildlife management domains. Besides these goals, different objectives have been fulfilled by discovering motifs recently. Hence, the work presented in Tanaka et al. (2005) proposed an algorithm devoted to discover motifs based on the minimum description length principle. The approach also allowed to obtain motifs from multi-dimensional time series data by using principal component analysis. Mueen et al. (2009) proposed in 2009 an exact algorithm to find time series motifs much faster than brute-force searching strategies do. Also, an approach based on tree-construction search to discover motifs in multivariate time series was proposed in Wang et al. (2010) and applied to sensory datasets. Since Stormo (2000) first reviewed strategies to find DNA motifs (meaningful base sequence patterns that identify binding sites responsible for transcription factors) in 2000, a large amount of algorithms have been developed. Thus, an ensemble algorithm attempting to discover regulatory motifs in DNA sequences was proposed in Hu and Kihara (2006). Another algorithm was proposed in Wijaya et al. (2007) which, given a set of sequences, executes mdifferent motif finders, each of them reporting nmotifs. Finally, Sharov and Minoru (2009) presented CisFinder, a software that generates a comprehensive list of motifs enriched in a set of DNA sequences and describes them with position frequency matrices. 2.3. Robust statistical methods to detect outliers The problem of a posteriori outliers detection in time series has been widely studied in the literature, and faced from many different points of view. In fact, the existence of even few outliers usually leads to inaccurate models and not satisfactory forecasts (Galeano et al., 2006), since they may deeply influence the estimates that classical methods propose (Carnero et al., 2007). For this reason, there is a large family of robust statistical methods (Rousseeuw and Hubert, 2011) that deal with outliers and, particularly, propose approaches to detect their existence in the datasets subjected to analysis. Gelper et al. proposed an adapted version of the classical exponential and Holt-Winters smoothing methodologies, providing them with robustness (Gelper et al., 2010). Another version of a robust multivariate exponential smoothing applied to time series can be found in Croux et al. (2010). Following with classical methods, a work that enhanced ARMA by adding robustness can be found in Muler et al. (2009), in which the authors succeeded in limiting the effect of outlying data to the time stamp in which they happen. Support vector machines (SVM) have also been adapted to deal with outliers. Actually an approach to model time series using robust SVM was proposed in Camps-Valls et al. (2004). In particular, the authors claimed that their proposal provides stable models and allows the analysis of models’ memory depth. Recently, a twosteps methodology was proposed in Chuang and Lee (2011) that combinedtheuseofrobustSVMtoremoveanomalousobservations, and non-robust SVM to obtain estimates from that reduced dataset. Many of these proposals have been implemented and freely distributed in software packages. But from all of them there are two that highlight. LIBRA (Verboven and Hubert, 2005) is a Matlab library for robust analysis, that contains (among others) robust covariance estimation, regression, principal component analysis, principal component regression or partial least squares, as well as methodologies to detect outlying observations in datasets. The TOMCAT toolbox Daszykowski et al. (2007), also developed in the Matlab environment, includes almost the same methods that LIBRA does, but also includes a graphical interface. 3. Fundamentals This Section first defines some terms in order to prevent possible misinterpretations in sensitive terms. Since the proposed methodology is based on an existing algorithm, this Section also provides a brief summary of the mathematical fundamentals underlying the PSF algorithm. Note that a more detailed explanation can be found in Martínez-Álvarez et al. (in press). 3.1. Definitions This work uses certain concepts –such as outlier or motif– that can be interpreted in many different senses, depending on the application or even the author. For this reason, this Section provides a formal definition for these sensitive terms. Definition 1 (Hourly time series). An hourly time series Tis a set of real-valued data in successive order, occurring every hour. In this work, T=[t 1 ,...,t p ], where pis the length of the time series and usually a multiple of 24. Definition 2 (Daily time series). From an hourly time series, a daily time series Dis formed by tuples in R 24 ,D=[d 1 ,...,d p/24 ], where d i =[t 24(i1)+1 ,...,t 24i ]. Definition 3 (Label). In this work, the term label is used to identify a set of possible categorical values. Thus, L¼fl 1 ;...;l K g, where Kis a pre-fixed number. Definition 4 (Sequence). A sequence Sis a set of labels occurring in successive order. In this work, S=[s 1 ,...,s q ], where qis the length of the sequence and s i 2L. Definition 5 (Outlier). Given a test set, an outlier is an observation which appears to be inconsistent with the rest of the data, relative to an assumed model (Everitt, 2006). Outliers are usually represented by a binary random variable v i , for i=1,...,pthat models their occurrence ( v i = 1 if it occurs, and v i = 0 otherwise), and by another real random variable z i , for i=1,...,pthat models their magnitude (Maronna et al., 2007). Although different outlier types can be found in the literature, only the additive outlier model is considered in this work, due to the nature of the studied data: t i ¼x i þ v i z i ;ð1Þ where t i the observed value, and x i the i th cleaned data modeled by any approach. Definition 6 (Motif). A motif M W is a sequence of Wconsecutive labels considered to occur just before an outlier, where Wis the pre-fixed length of the sequence. In addition, M W is a subsequence found in Sand, consequently: M W ¼½s 0 1 ;...;s 0 W , where s 0 i 2L. 3.2. Time series forecasting: the PSF algorithm The PSF algorithm is a general-purpose time series forecasting algorithm whose main feature is that it only makes use of certain labels –obtained by means of a clustering process– to forecast arbitrary horizons of prediction. However, the output is not composed by labels but by real values. PSF can deal with an arbitrary number of samples per day. However, specifically in this paper, the time series considered consists of twenty-four samples per day. That is, given the hourly values up to day ifor a time series, the PSF algorithm provides the 24 hourly values corresponding to day i+ 1. Formally, let d i 2R 24 be a vector that comprises the 24 hourly values of a certain day, i. First, PSF applies clustering techniques to such data in order to assign a label to each day. Formally, it uses a function F K that assigns a label l i 2Lto the values d i 2Dof each day by means of a clustering process, F K :D!L, that is, every 24 h are identified by a label. Once Kis fixed, this process transforms the daily time series Dinto a sequence of labels S, thus discretizing the original data. Let l i be the label assigned to the day iobtained by means of the application of a clustering technique. Let S i W be the labels’ subsequence of Wconsecutive days, from day ibackward: S i W ¼½l ið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. LetW ⁄ be the length of the window determined by PSF. For a day i and length of window W ⁄ , the PSF algorithm searches for the subsequences of labels which are exactly equals to S i W  in the dataset, providing the equal subsequences set, ES, defined by the equation, ESði;W  Þ¼ days j2Dsuch that S j W  ¼S i W  no ð3Þ It is worth remarking that if no subsequence equal to S i W  was found in the dataset, that is, ES(i,W ⁄ )=;, the length of the window would decrease by one unit, W 0 =W ⁄ 1, and the PSF would search for subsequences equal to S i W 0 . This process may be repeated until any subsequence is found, that is, ES(i,W)–;. Therefore, the W ⁄ consecutive labels that precede the day to be predicted are extracted and searched for in the historical data. Once all occurrences of S i W  are found, the 24 hourly values of the day i+ 1 are predicted by averaging the real values found immediately after each S i W  match. Mathematically, ^ d iþ1 ðW  Þ¼ 1 #ESði;W  ÞX j2ESði;W  Þ d jþ1 ð4Þ Finally, the daily error for any day iis defined by: e day ði;W  Þ¼j ^ d i ðW  Þd i jð5Þ 4. Outlier forecasting in time series This section explains the methodology proposed to improve the forecasting process provided by the PSF algorithm. The discovery of motifs is included in the aforementioned algorithm as a crucial step for forecasting the occurrence of outliers and, then, providing accurate estimates for such anomalous observations. The value of two parameters had to be determined in the PSF process: The number of clusters Kand the length of the window W. With regard to K, the new approach acts exactly the same as what was proposed in the original PSF, that is, it applies three well-known validity indices –Silhouette, Dunn and Davies–Bouldin– and determines the optimal number of clusters by means of a majority vote system. On the other hand, the nfold cross-validation is used to obtain the optimal value of W. Twelve folds have been created in this work (n= 12) for all the datasets, where each fold represents a month. Therefore the training set consists of one year. The 12 fold cross–validation is then evaluated. The forecasting errors are calculated in every fold by varying the length of W. For each window size W, the monthly errors are denoted by e month (W) and are calculated as follows: e month ðWÞ¼ 1 #month X i2month e day ði;WÞð6Þ for W=1,...,W max and W max = 10, since no longer sequences were found in daily time series. Then, the average errors are calculated for each window size as follows,  eðWÞ¼1 nX month e month ðWÞð7Þ where n= 12 and month ={Jan,...,Dec}. The W ⁄ selected is the one that minimizes the average error corresponding to the 12 folds (months) evaluated. W  ¼argminf eðWÞg with W¼1;...;W max ð8Þ It is now –just after the training step and before the prediction process– that the discovery of motifs plays a crucial role, as it attempts at forecasting the occurrence of anomalous days in time series. Specifically, the motifs to be found are those which generate a prediction error greater than the average error in the cross–validation process. Therefore, a set of days ibelonging to the training set (TS) that satisfies e day ði;W  Þ> eðW  Þis constructed. This set, CS or candidates set, gathers all the candidate days to be preceded by a sequence that will eventually be a motif. Formally, CS ¼fi2TS such that e day ði;W  Þ> eðW  Þg ð9Þ Nevertheless, not all the sequences that precede these candidates have the same probability to be eventually considered as outlier precursors, since the associated errors range from values close to the mean error (these candidate sequences should be eventually discarded by the approach) to significantly high values. For this reason, each candidate is co-labeled by using clustering techniques, more specifically, the K-means algorithm. The decision on how many clusters have to be created is always an open question and many indices could be used. However, it is worthless to have a large number of clusters and therefore, only three conceptual classes will be created: A class that gathers days with low errors (C l ) or nearest errors to the  eðW  Þ, a class containing medium errors (C m ) and, finally, a class devoted to identify high errors (C h ) or farthest errors to the  eðW  Þ. Thus, for all days i2C l ,j2C m and k2C h the following inequalities are fulfilled: e day ðk;W  Þ>e day ðj;W  Þ>e day ði;W  Þ> eðW  Þð10Þ Fig. 1 illustrates an imaginary error distribution in the TS, determining the candidates that will eventually form CS as the union of candidates in C l ,C m and C h . In other words, CS =C l SC m SC h . A priori reasoning reveals that those candidates belonging to C h must be more probable to be preceded by motifs preceding outliers than the candidates in C l or C m . Results corresponding to each cluster of data will be separately analyzed in Section 5. The next step consists of computing the sequences of labels occurring before the candidates in order to determine which sequences will be considered outlier precursors, since not all these sequences will be motifs. Hence, the approach has to decide whether the sequences preceding the candidates are motifs responsible for outliers or not. In particular, the sequences that only appear before the candidates will be considered motifs preceding an outlier occurrence. That is, if a sequence of labels preceding a candidate is found prior to any other day that does not belong to the CS, the sequence is discarded and not considered to be a motif. Thus, a set of motifs MS is defined by: MS ¼S i W  such that i2CS and ESði;WÞ#CS no ð11Þ Note that MS can be also expressed as follows: MS =M l SM m SM h , where M l ,M m and M h are the discovered motifs associated to the sequences found in the classes C l ,C m , and C h respectively. Fig. 4 depicts the process of discovering the sequences of labels that represent the motifs preceding outliers. This figure shows an illustrative example in which six clusters were created, K=6 (labels are digits 1 to 6). Therefore, the time series appears discretized, making only use of six different labels to be assigned one to each day. The labels in bold refer to those days initially included in CS, that is, those that obtained forecasting error greater than the average. By contrast, the labels followed by a bullet are those days that do not belong to CS but are preceded by a sequence equal to another that precedes a candidate. Then, the three labels preceding each day in CS (W ⁄ = 3 in this example) are extracted. Finally, the cases in which the sequences before the candidates also appeared before any day not belonging to CS were discarded. Otherwise, these sequences are considered motifs or pattern sequences preceding an especially unexpected value. For this particular example, note that only the sequence {1,5,4} would have been considered as motif. The goals are now: (i) to forecast the occurrence of an outlier and (ii) provide accurate estimation for it. Thus, once MS is 50 100 150 200 250 300 350 0 5 10 15 Days Error (%) Error distribution in TS Cl Ch Cm MEAN Fig. 1. Illustrative distribution of candidates in C l ,C m and C h . constructed, the general scheme of forecasting is as follows. Original PSF just extracted S W  and searched for it in the historical data. But now, before its search, it has to be determined if this pattern sequence matches any of the motifs forming MS. Given this situation, two cases may arise: (1) S W  does not match any of the motifs in MS. The approach determines that the day to be forecasted is not an outlier and it would continue with normal PSF procedure. That is, ^ d iþ1 ðW  Þis calculated as PSF does from Eq. (4). (2) S W  matches any of the motifs in MS. The approach determines that the day to be forecasted is an outlier. In this case the challenge is to provide accurate estimations for these observations, i.e. to provide accurate values for z i . To fulfill this goal, a simple strategy is proposed: To average the value of the outliers found by the proposed approach in the historical data. Formally, the set of outlying days from the training set is defined by: OS ¼fi2TS such that S i W  2MSgð12Þ Then, the forecast for the outlier is based on the a posteriori detected outliers in the training set: ^ d iþ1 ðW  Þ¼ 1 #OS X i2OS d iþ1 ð13Þ where #OS is the number of outliers detected by the approach in the historical data (the number of elements in OS), and ^ d i the values of the time series for these outlying days that form OS. Fig. 2 illustrates the entire process of prediction when the discovery of motifs is included in the PSF algorithm. Note that this step has to be performed immediately after the clustering (creation of the sequence of labels) and before the forecasting. In addition, the steps corresponding to discovery of motifs and prediction are further detailed in Fig. 3. Finally, the pseudocode of the proposed methodology is presented in Fig. 5, and that of the discovery of motifs process in Fig. 6. 5. Results This section presents the results obtained by the application of the proposed methodology to six time series. The motifs discovery process to six real-world time series is described in Section 5.1. Then, a statistical analysis has been carried out to determine the validity of the assumptions made when forecasting the occurrence of outliers in the time series. This analysis can be found in Section 5.2. Finally, to compare the results obtained, Section 5.3 reports average forecasting errors of the new methodology and other techniques. 5.1. Motifs discovery in real-world time series The discovery of motifs on real-world time series is now described. In particular, six public electricity-related (three of prices and three of demand) time series have been considered to show that the proposed methodology properly works on different datasets. Thus, the new approach has been applied to the Spanish (OMEL), New Yorker (NYISO) and Australian (ANEM) markets, whose data are available on-line in Spanish Electricity Price Market Operator (http://www.omel.es), The New York Independent System Operator (http://www.nyiso.com) and Australia’s National Electricity Market (http://www.nemmco.com.au), respectively. The forecasting process is applied to the year 2006 for the three markets, with a historical data of one year and with a horizon of prediction of one month. As twelve months are going to be evaluated for each market, the methodology is going to be tested on 72 datasets. Given this situation, every time a month is forecasted the training set changes. For instance, when January 2006 is forecasted, the training set comprises the whole year of 2005. However, when February 2006 is forecasted the historical data ranges from February 2005 to January 2006, and so on. These changes in the training set involve changes in the configuration of PSF. First of all, both Kand Whave to be determined according to the methodology presented in Section 4.Table 1 summarizes the values of these parameters for the six markets in the year 2006. The results after the motifs extraction step for the three markets are summarized inTables2–4.Thatis, a summaryofall encountered classes, sequences and motifs can be found in theses Tables for the Spanish,NewYorkerandAustralianmarkets,respectively.However, onlyelectricity prices resultsof January 2006 for theAustralianmarket are now described as the explanation for the remaining eleven months for each year and market is similar. Therefore, all the comments about the results provided below refer to prices shown in Table 4. First, the parameters to be set in the PSF are equal to: (K,W) (3,6), according to Table 1. The CS can be now constructed. For thispurpose, the  eðWÞ(seeEq.(7))has to beconsidered since the candidates are those days belonging to the training set (January to December 2005) that obtained an error greater than  eðWÞ.Thevalue of the mean error, calculated according to the methodology in Section 4is  eð6Þ¼5:81%. Therefore, CS would be formed by all days in 2005 with forecasting error greater than 5.81%. Now the three classes are constructed by applying K-means, with K= 3 as mentioned in Section 4. The three clusters are defined as: C l is the class that contains the candidates days with error from 5.81% to 7.13%, C m the one that gathers the candidates with errors ranging from 7.13% to 9.47% and C h the class that contains the candidates with errors greater than 9.47%. The error distribution, according to these three clusters, is shown in Fig. 7. Fromthe365 daysof2005thatcomprisethetrainingset,137(see Table 4, row 1: #C l +#C m +#C h = 101 + 32 + 4 = 137) had an error greater than 5.81% so the constructed CS contains 137 candidates days. Once the candidates are selected, the number of different sequences that generated them are considered. From the candidates in C l , 5 different sequences were found (S l = 5); from the candidates in C m ,3(S m = 3) and from the candidates in C h ,2(S h = 2). This fact involves that from all the K W = 729 possible sequences, only 10 caused errors greater than the average. Note that there were situations in which a particular sequence appeared before different candidates that belong to different classes. For these cases, only the sequence that appeared in the class with higher associated error (C l was devoted to include days with lower errors, C m to medium errors and C h to higher ones) was counted. Finally, the number of motifs that identify outliers are determined. From the sequences C l , only one appeared exclusively as Data Clustering MOTIFS DISCOVERY Prediction Labeled data Motifs Insert predicted sample More days? Yes End No Fig. 2. Illustration of the proposed methodology. an outlier precursor, M l = 1. With reference to the sequences in C m , one out of three, M m = 1. Last, both sequences in C h were exclusive, M h =2. Fig. 8 is provided to determine the usefulness of dividing the CS into three groups. These histograms show the motifs distribution along with C l ,C m and C h and, to be precise, the relation between different sequences and motifs for each class and market, expressed as %. Thus, each bar is calculated by dividing the number of different sequences that precede the candidates and the number of motifs that are finally selected. As it is possible to observe, the probability that a sequence becomes a motif is directly related with the error associated to the candidate to which it precedes. For instance, the percentage of motifs for the ANEM’s electricity price time series are 20.26% for the sequences in C l , 29.82% for the sequences in C m and 52.94% for the sequences in C h . The motifs found in January 2006, represented as a numerical sequence of labels, are shown in Tables 5 and 6 corresponding to prices and demand, respectively. For instance, note that the first motif found in ANEM’s prices is M 1 l ¼f1;3;2;2;3;1g. Each label is represented by one of the Kclusters generated during the training of the PSF (where K= 3 in this case) and identifies 24 h. As the length of the window was set to W= 6, these six labels in fact represent 144 h. Figs. 9 and 10 illustrate the most representative motifs found in all the markets when forecasting January 2006. Actually, these motifs represent the time series values that will precede an outlier. As each motif was represented by Wconsecutive labels (see Definition 6), these figures depict the values associated to every label, which have been obtained by means of clustering techniques. Construct CS Clustering CS Find motifs Candidates DISCOVERY OF MOTIFS {Cl,C m,C h} Search for S in historical data PREDICTION Motifs set, MS Window's length, W* any motif? NO YES Forecasts S matches w* w* Search for motifs Forecast data from Eq. (13) in historical data Forecast data from Eq. (4) Forecasts (Outlier occurrence) Motifs set, MS Window's length, W * Training set, TS Number of clusters, K Set W W* Fig. 3. Detail of discovery of motifs and prediction steps. Fig. 4. Illustrative example of motifs discovery. Fig. 5. A general scheme of the proposed methodology. Regarding the electricity prices, note that for the Spanish market, six motifs were found; six for New York, and four for the Australian market. As for the electricity demand, the number of motifs found were seven, five and four for the Spanish, New York and Australian markets, respectively. In actual fact, the curves in these figures represent the average evolution of the five (OMEL), four (NYISO) and six (ANEM) days prior to an outlier forecast in prices time series and the average evolution of the three (OMEL and NYISO) and four (ANEM) days prior to an outlier forecast in demand time series. Fig. 6. Pseudocode of the motifs discovery. Table 1 Setting the PSF. OMEL refers to the Spanish market, NYISO to the New Yorker market and ANEM corresponds to the Australian market. Month Prices Demand OMEL NYISO ANEM OMEL NYISO ANEM KWKWKWKWKWKW January 4 5 5 4 3 6 7 3 4 3 6 4 February 4 5 5 3 3 6 7 4 3 5 6 5 March 4 5 5 4 3 6 7 3 4 4 4 3 April 4 5 5 4 4 6 6 4 4 5 4 4 May 45 63 46 54 44 64 June 6 4 5 3 3 6 6 3 4 5 5 5 July 5 5 6 4 3 5 6 4 3 5 6 4 August 6 4 6 3 4 6 5 4 4 3 5 4 September 6 4 5 3 3 6 7 3 5 4 6 5 October 6 4 5 4 3 6 5 4 5 4 5 4 November 6 4 5 4 3 6 5 3 4 4 6 4 December 5 5 5 3 3 5 6 3 5 3 4 6 Table 2 Motifs distribution for Spanish markets. Month Prices Demand #C l (S l )[M l ]#C m (S m )[M m ]#C h (S h )[M h ]#C l (S l )[M l ]#C m (S m )[M m ]#C h (S h )[M h ] January 98(7)[1] 25(4)[3] 8(2)[2] 103(6)[2] 35(6)[2] 6(3)[3] February 87(6)[0] 31(7)[2] 5(2)[2] 121(8)[1] 22(5)[3] 3(2)[1] March 73(5)[2] 16(3)[1] 8(1)[1] 99(7)[1] 30(7)[4] 6(1)[0] April 103(9)[1] 30(6)[1] 6(3)[2] 83(4)[0] 17(5)[2] 5(3)[3] May 65(4)[0] 51(6)[0] 10(4)[2] 101(6)[2] 26(8)[3] 5(3)[3] June 97(6)[0] 38(5)[2] 4(0)[0] 124(7)[3] 27(5)[1] 7(1)[1] July 180(8)[2] 27(3)[1] 12(5)[3] 97(8)[2] 13(3)[0] 8(4)[3] August 101(8)[3] 26(5)[0] 9(4)[2] 113(11)[4] 29(6)[3] 11(3)[2] September 110(8)[1] 25(5)[3] 5(0)[0] 105(9)[2] 19(5)[4] 3(2)[2] October 108(7)[1] 23(4)[2] 6(1)[0] 98(8)[3] 32(3)[2] 4(1)[1] November 120(9)[1] 40(6)[0] 6(3)[1] 99(9)[2] 20(5)[3] 9(3)[3] December 169(10)[2] 38(9)[3] 10(3)[1] 114(8)[2] 24(3)[2] 12(2)[0] Fig. 11 shows ANEM’s electricity prices for July 2006. This month is illustrated since it presents large outliers. As for this month W= 5 (see Table 1), the number of labels considered for outlier occurrence forecasting is five. Also, the number of motifs found in TS were five (see Table 4): M 1 l ¼f1;2;3;3;1g;M 2 l ¼f2;2; 3;1;3g;M 3 l ¼f3;1;1;1;1g;M 4 h ¼f2;1;3;3;1gand M 5 h ¼f2;1;1; 1;3gThus, grey bars represent outliers a posteriori detected by means of the robust statistical method presented in Gelper et al. (2010), that is, days 3, 13, 18 and 22 July were the outliers identified. It can be appreciated that sequences M 3 l ;M 4 h and M 1 l are three of the motifs found in TS that eventually preceded days 13, 18 and 22 July, respectively. Furthermore, 3rd July was preceded by the motif M 5 h . Only the last two labels of this motif (1,3) are depicted in Fig. 11 because the first three labels (2,1,1) correspond to days in June. Finally, note that one of the motifs found in TS, M 2 l , did not occur when forecasting July 2006. 5.2. Evaluation of outlier occurrence forecasting Once all the MS have been constructed, the proposed method predict a priori if the day to be forecasted will be an outlier. This Section is devoted to statistically quantify the outlier occurrences Table 3 Motifs distribution for New York markets. Month Prices Demand #C l (S l )[M l ]#C m (S m )[M m ]#C h (S h )[M h ]#C l (S l )[M l ]#C m (S m )[M m ]#C h (S h )[M h ] January 101(8)[3] 34(3)[1] 12(2)[2] 96(9)[1] 28(3)[1] 14(4)[3] February 92(11)[2] 36(4)[1] 14(3)[2] 88(9)[3] 14(2)[1] 7(2)[2] March 89(7)[2] 45(5)[1] 11(2)[1] 101(9)[2] 33(2)[0] 9(3)[2] April 110(13)[3] 21(4)[1] 9(5)[4] 93(8)[2] 29(3)[2] 11(4)[2] May 121(7)[2] 31(2)[1] 6(3)[1] 114(9)[3] 30(4)[2] 13(6)[3] June 142(5)[1] 32(0)[0] 7(0)[0] 103(7)[0] 26(3)[1] 10(4)[2] July 92(10)[3] 41(5)[1] 18(4)[4] 86(5)[1] 29(4)[2] 8(3)[2] August 84(7)[2] 39(6)[2] 9(6)[5] 76(7)[4] 14(2)[2] 4(1)[1] September 107(7)[0] 40(4)[0] 10(1)[1] 84(8)[1] 21(4)[1] 8(3)[3] October 141(9)[2] 28(3)[0] 4(0)[0] 95(6)[0] 32(4)[2] 9(3)[1] November 99(12)[3] 32(8)[2] 15(4)[3] 115(8)[3] 38(6)[2] 15(6)[2] December 87(8)[1] 44(3)[1] 8(0)[0] 109(8)[3] 29(4)[2] 12(4)[3] Table 4 Motifs distribution for Australian markets. Month Prices Demand #C l (S l )[M l ]#C m (S m )[M m ]#C h (S h )[M h ]#C l (S l )[M l ]#C m (S m )[M m ]#C h (S h )[M h ] January 101(5)[1] 32(3)[1] 4(2)[2] 131(7)[0] 46(7)[2] 5(3)[2] February 165(5)[0] 25(6)[1] 9(1)[1] 115(6)[1] 56(7)[3] 6(3)[3] March 133(8)[2] 13(2)[0] 3(0)[0] 92(6)[0] 72(9)[1] 7(3)[3] April 190(13)[3] 8(2)[0] 11(4)[2] 102(7)[2] 64(5)[0] 5(4)[3] May 187(17)[5] 13(4)[2] 6(3)[1] 123(7)[2] 41(3)[1] 11(4)[0] June 169(12)[5] 22(5)[1] 3(3)[1] 183(16)[3] 33(3)[0] 4(4)[2] July 172(22)[3] 9(0)[0] 2(2)[2] 117(9)[1] 40(5)[0] 8(5)[3] August 142(20)[2] 34(3)[2] 6(4)[2] 139(11)[2] 41(4)[3] 7(3)[1] September 102(18)[4] 20(8)[3] 12(2)[1] 110(8)[2] 27(3)[2] 9(3)[3] October 81(13)[2] 53(13)[4] 9(2)[0] 142(12)[3] 42(8)[3] 5(1)[1] November 112(9)[1] 43(7)[1] 14(5)[4] 137(12)[2] 38(7)[3] 7(3)[2] December 121(11)[3] 39(4)[2] 13(6)[2] 134(10)[4] 29(5)[2] 8(4)[2] 0 50 100 150 200 250 300 350 0 2 4 6 8 10 12 14 16 18 20 Days Error (%) 9.47% 7.13% 5.81% Ch Cm Cl Fig. 7. Forecasting errors in TS distributed in C l ,C m and C h obtained by applying the K-means to the candidate set CS with K=3. F. Martínez–Álvarez et al. /Pattern Recognition Letters 32 (2011) 1652–1665 1659 properly forecasted by the proposed methodology. Thus, the quality parameters used to evaluate the accuracy of the approach are first introduced in Section 5.2.1 and, then, the conducted statistical analysis is reported in Section 5.2.2. 5.2.1. Parameters of quality The parameters used to assess the accuracy of the approach are now introduced. A posteriori analysis has been carried out to determine the existence of outliers in all examined series. In particular, the robust method proposed in Gelper et al. (2010) (hereafter called RHW for simplicity) to detect outliers in time series has been considered. Hence, a forecast of an outlier occurrence is said to be properly made by the proposed approach, if RHW also points the observation as anomalous. Thus, a priori forecasting (the proposed approach in this work) is compared to a posteriori detection (the method proposed in Gelper et al. (2010)). Note that the authors in Gelper et al. (2010) determined that outliers are those data that do not fulfil any of the two bounds they define: UB t ¼^ y t þ2^ r ð14Þ LB t ¼^ y t 2^ r ð15Þ where UB t and LB t are the upper and lower bounds respectively, ^ y t is the fitted value, and ^ r is the standard deviation of the regression residuals they obtain. Hence, in subsequent equations, true positives or TP is the number outlier occurrences properly forecasted, that is, the number of days preceded by a motif in MS that are outliers according to RHW; true negatives or TN is the number of days that were not preceded by a motif in MS and were not considered outlier by RHW either; false positives or FP is the number of days preceded by a motif in MS that were not considered outlier by RHW; and false negatives or FN is the number days not considered outliers (not preceded by any motif in MS) and eventually considered outliers by RHW. According to these definitions, the sensitivity is the probability that a motif discovered precedes a real outliers. Its formula is defined as follows: Sensiti v ity ¼TP TP þFN ð16Þ Another relevant parameter is the specificity, or the ratio of sequences preceding the day to be forecasted properly discarded by the approach. The mathematical expression is: Specificity ¼TN TN þFP ð17Þ The positive predictive value (PPV) is the probability that a forecasted outlier is indeed a real one. Its formula is: PPV ¼TP TP þFP ð18Þ 0.00% 10.00% 20.00% 30.00% 40.00% 50.00% 60.00% 70.00% 80.00% 90.00% OMEL ANEM NYISO CmCh Cl 0.00% 10.00% 20.00% 30.00% 40.00% 50.00% 60.00% 70.00% 80.00% 90.00% OMEL ANEM NYISO ClCmCh Fig. 8. Motifs distribution in C l ,C m and C h . Table 5 Motifs found in TS for the electricity price when forecasting January 2006. Market Motif OMEL M 1 l ¼f1;1;3;4;4g M 2 m ¼f3;4;4;1;3g M 3 m ¼f3;4;4;2;1g M 4 m ¼f2;3;4;4;1g M 5 h ¼f1;1;1;2;3g M 6 h ¼f3;2;3;3;4g NYISO M 1 l ¼f1;3;5;4g M 2 l ¼f2;2;3;1g M 3 l ¼f5;1;1;3g M 4 m ¼f3;4;2;5g M 5 h ¼f3;2;4;1g M 6 h ¼f3;1;5;1g ANEM M 1 l ¼f1;3;2;2;3;1g M 2 m ¼f2;1;3;3;3;1g M 3 h ¼f3;2;1;1;2;3g M 4 h ¼f3;3;3;1;1;2g Table 6 Motifs found in TS for the electricity demand when forecasting January 2006. Market Motif OMEL M 1 l ¼f6;4;4g M 2 l ¼f7;1;7g M 3 m ¼f5;5;4g M 4 m ¼f4;7;5g M 5 h ¼f4;4;7g M 6 h ¼f6;5;6g M 7 h ¼f7;7;2g NYISO M 1 l ¼f2;3;3g M 2 m ¼f4;3;4g M 3 h ¼f1;3;2g M 4 h ¼f3;4;1g M 5 h ¼f2;2;3g ANEM M 1 m ¼f6;5;3;4g M 2 m ¼f2;6;3;4g M 3 h ¼f2;3;1;4g M 4 h ¼f4;6;6;5g