UniFIeD Univariate Frequency-based Imputation for Time Series Data

Schriftenreihe CIplus, Band 5/2013 Herausgeber: T. Bartz-Beielstein, W. Konen, H. Stenzel, B. Naujoks UniFIeD Univariate Frequency-based Imputation for Time Series Data Martina Friese, J¨org Stork, Ricardo Ramos Guerra, Thomas Bartz-Beielstein, Soham Thaker, Oliver Flasch, Martin Zaefferer UniFIeD Univariate Frequency-based Imputation for Time Series Data Martina Friese, Jorg Stork, Ricardo Ramos Guerra, Thomas Bartz-Beielstein, Soham Thaker, Oliver Flasch, Martin Zaefferer Faculty for Computer and Engineering Sciences Cologne University of Applied Sciences, 51643 Gummersbach, Germany [email protected] Abstract. This paper introduces UniFIeD, a new data preprocessing method for time series. UniFIeD can cope with large intervals of missing data. A scalable test function generator, which allows the simulation of time series with different gap sizes, is presented additionally. An experimental study demonstrates that (i) UniFIeD shows a significant better performance than simple imputation methods and (ii) UniFIeD is able to handle situations, where advanced imputation methods fail. The results are independent from the underlying error measurements. 1 Introduction Missing data is a well-known problem in nearly every real-world time series. Sensors may fail, data might get lost during transfer, or measurements are simply missing. Although this problem is well-known, many standard time-series prediction and analysis methods rely on complete data. During the last decades, several ways have been developed to tackle missing data. Several of these methods are applicable to small gap sizes only. A common suggestion, which is available in several software packages, is the imputation of mean values. This approach can destroy inherent data structures and may worsen the statistical modeling, resulting in large prediction errors [4]. Imputing missing elements from regression or analysis of variance (ANOVA) are usually better. More advanced methods, also from computational intelligence, for data imputation were developed in the context of univariate (linear, spline, and nearest neighbor interpolation), multivariate (regression-based imputation, nearest neighbor, self-organizing map, multi-layer perceptron), and hybrid methods of the previous by using simulated missing data patterns. A small study, which discussed the applicability of these methods to air quality data sets, was performed by Junninnen et al. [7]. Single imputation methods, i.e., filling in precisely one value for each missing one, can be distinguished from multiple imputation methods. The latter generate multiple simulated values for each missing value. Our study was initialized by an real-world task and focuses on the applicability of univariate methods. It was motivated by a real-world problem, because we received time-series data with large gaps from one of our industrial partners. These data should be used for time-series predictions, where the methods 2 Friese et al. of choice require complete data. Since simple imputation methods failed completely in our setting and the advanced methods did not show the expected performance, we decided to develop a new imputation method. The new method, entitled univariate frequency-based imputation for time series data (UniFIeD) outperformed the existing methods. The success in field settings right from the start motivated a first analysis and gave reason for performing an experimental study. Focussing on methods for large intervals of missing data, seasonal data, especially time series data, we consider the following scientific goals: (G-1) Which method generates the smallest imputation error? (G-2 ) What is the influence of data pre-processing methods on the performance of forecast methods? Based on these goals, we are interested in developing an automated and robust procedure for data pre-processing, which can be implemented easily. To generate scientifically significant results, we will proceed as follows. First we generate instances based on the real-world data. Then we run the imputation methods. Next, the errors based on different error measurements are determined. Finally, their prediction errors are reported and compared on different error measurements. As a future step to increase the plausibility of our findings, we are planning to perform predictions with different state-of-the-art methods. This paper is structured as follows: First, the real-world data is described in Sec. 2. Pre-processing methods and the univariate frequency-based imputation for time series data (UniFIeD) are introduced in Sec. 3. The prediction models, which were used in the final comparison, are described in Sec. 4. Error measures, which play a crucial role in our study, are presented in Sec. 5, the two different experiments are introduced in Sec. 6. Results are presented in Sec. 7, our findings are discussed in Sec. 8. The paper concludes with a summary and an outlook in Sec. 9. An R version of the program code used in this study, is freely available for download and will be compiled as an R package [8]. 2 Data 2.1 Missing Data We consider three types of data: y∗denotes the underlying (latent and complete) data, yis the observed data, and ˆyis the imputed data. To evaluate the performance of an imputation method, criteria have to be defined. The imputation performance depends (at least) on two characteristics: (a) the structure of missing data pattern and (b) the amount of missing data. If the probability of missing data does not depend upon the observed or the unobserved data, then these data are called missing completely at random (MCAR)[9]. There is no predictive power in the observed values y, if the missing value process is MCAR. In general, the structure of missing data in our projects is MCAR. The simulation of missing data pattern randomly will be described in Sec. 2.3. Univariate Frequency-based Imputation for Time Series Data 3 2.2 The Datasets Our study is based on real-world data. The experiments are based on energy consumption time-series data supplied by GreenPocket GmbH. The data was recorded by two independent smart metering devices, installed at a local commercial customer. Some data points are missing due to measurement or transmission issues, which is a common situation in real-world settings. The data provided by GreenPocket GmbHis a series of timestamp and meter reading pairs taken quarter hourly. Timestamps are given an ISO 8601 derived date/time format, meter readings are given in kilowatt hours (kWh). The energy consumption time series data was recorded at the same time interval by two independent smart metering devices resulting in the two data sets series_meter1 and series_meter2. Both energy consumption time series datasets contain a total of 8548 entries starting at 2010-12-06 23:15:00 and ending at 2011-03-06 00:00:00 which makes a total time interval of more than 12 Weeks. The complete time series data set series_meter1 is shown in Figure 1, whereas Figure 2 shows only the last two weeks of the same data set. [1] "English_United States.1252" [1] "German_Germany.1252" 0 2 4 6 8 Dec 15 Jan 01 Jan 15 Feb 01 Feb 15 Mar 01 energy consumption (kWh) Fig. 1. plot of all the observations given in series_meter1 Visual inspection of the data reveals daily periods, while weekly periods are detectable, but not as clearly defined. Also time intervals with missing data can be clearly seen. Having a closer look at the missing data values, contained in both time series, reveals that there are altogether twelve gaps. Mostly smaller gaps of length one, but also larger gaps up to the size of 385 missing observations. 2.3 Test Instance Generation Since missing data already occurs in the real-world test data, we are able to determine a realistic distribution of the gap sizes and frequencies. This includes 4 Friese et al. [1] "English_United States.1252" [1] "German_Germany.1252" 0 2 4 6 Feb 21 Feb 23 Feb 25 Feb 27 Mar 01 Mar 03 Mar 05 energy consumption (kWh) Fig. 2. excerpt of the observation from series_meter1 showing only the last two weeks of the data set the determination of two parameters: (a) distribution of the gap sizes gapsize and (b) the total amount of missing data, i.e., gappercentage. The gapsize distribution can be estimated from real-world data as follows. First, histogram plots were used for visual inspection of the gap sizes. Visual inspection suggests an exponential distribution of gapsize with smaller gaps appearing more frequently than larger gaps. The parameter λfrom the density of the exponential distribution f(t) = λexp(−λt) is estimated from the real-word data. Now we are able to generate random gaps with reasonable sizes that are in correspondence with real-world data. In a second step, we determine the percentage of missing data, gappercentage. Here, we consider values between 5 and 30. The generation of a single test instance than works as outlined in algorithm 1. 3 Pre-processing Methods 3.1 Existing Methods Existing imputation methods can be partitioned into two groups: the first group included basic methods [7], that do not use complex computations to determine the imputed values. A second group uses sophisticated techniques for imputations. We will present the basic methods first. Basic Imputation Methods Mean imputation is an often used method because of its simplicity. The missing data is replaced by the mean of the nonmissing observed data. ˆy=y, (1) where ˆyis the imputed value, ytare the observed values, and ydenotes the sample mean. Linear interpolation uses the start and end point of a gap to Univariate Frequency-based Imputation for Time Series Data 5 input : Time Series t input :gappercentage countData = number of observations ! = NA in t countDrop =countData ×gappercentage/100; repeat draw random gapsize from exponential distribution; until countDrop reached; countGaps = number of gaps drawn; remainingData =countData -countDrop; generate random partition of countGaps+1 parts summing up to the size of remainingData;(assuring that the generated time series neither starts nor ends with a gap and has at least one data point between two gaps.) t∗=t; for i= 1 →countGaps do position = sum(partitions[1:i] + sum(gapsize [1:i-1])); remove data from s∗from position to position+gapsize [i]; end output:t∗ Algorithm 1: Generation of a test instance construct a straight line. ˆy=yt1+k×(t−t1)with k =yt2−yt1 t2−t1 (2) yt1and yt2are the start and end values of the gap, while t1and t2are the start and end time values. x is the current time value. Nearest neighbors uses the start and end points of a gap as estimates for the imputation. ˆy=yt1if t < t1+t2−t1 2 ˆy=yt2if t > t1+t2−t1 2 (3) Advanced Imputation Methods The second group of imputation methods uses advanced techniques. The mice (Multivariate Imputation by Chained Equations) package specializes on multiple imputation methods [1]. The methods works best on multivariate data and no method applicable to the smart metering data was found to deliver good results. The zoo packages provides a methods na.structTS uses a generic function for filling NA values using seasonal Kalman filter [10]. Finally the Amelia II package can be mentioned here [3]. It was not able to find suitable values for the time interval from the Smart Metering data set. These packages obtain very good results, if multivariate time-series data were available. 3.2 Univariate Frequency-based Imputation The UniFIeD method proposed in this work relies on an automated estimation of time-series frequencies. 6 Friese et al. Estimating time-series frequencies automatically For the estimation of the frequencies contained in the data, the auto correlation function (acf) is used. The algorithm works as described in Algorithm 2. input : Time Series t determine acf values via auto correlation function on t; remember autocorrelation values from acf; remember related lags from acf; repeat reduce autocorrelation values to its peaks remember related lags; determine frequency from lags; via the frequency of the distances from one peak-lag to another until no new frequency found; output: all frequencies found Algorithm 2: Estimation of an underlying frequency using the auto correlation function. For better illustration, figure 3 shows the auto correlation function for series_meter1. Both, daily and weekly periods are clearly recognizable. The set of peaks that are considered as indicators for the data’s underlying frequency, are marked with small circles. In the second iteration of the algorithm, this set of peaks is reduced to a smaller set, which is marked with filled circles. After two iterations, the algorithm stops since there is no lower frequency in the data. How UniFIeD works UniFIeD is valuable for univariate time series data that presens a pattern or a seasonal effect. The algorithm takes advantage of this seasonal effect to find correlated patterns of the missing window to develop the imputation. Using the frequency estimated with Algorithm 2, we proceed to impute the missing values from the time series data set t. UniFIeD was developed not only to consider single missing points ˆybut full missing time windows of any size. The basic idea is to iteratively look for the next missing point ˆyof the time series tat the time moment tmand, using the frequency fand the number of similar windows kto search in, gather the amount of 2knon missing points ythat correlate to the time moment tmof the found missing point and form a vector ts. Once this vector tsis formed, and depending on the user’s request, the appropriate method to calculate the value to impute into the time series tis selected, where the available options to determine the new value are the mean, median, maximum, or minimum values of vector ts. Figure 4 shows an example of a missing value ˆy, the selected correlated values in the vector tsand the different methods used to impute the new value of a randomly selected missing window for the purposes of illustration. The UniFIeD algorithm works as follows: Univariate Frequency-based Imputation for Time Series Data 7 0 1000 2000 3000 4000 5000 6000 −0.4 −0.2 0.0 0.2 0.4 0.6 0.8 1.0 Lag ACF Fig. 3. Plot of the auto correlation function output on meter_series1. It shows the auto correlation value for each lag. Values, considered as peaks indicating the frequency, are marked with small circles. Note that the final step of Algorithm 3 is to indicate whether any NA were imputed into time series t. If this is fact, the algorithm will look for the left neighbor and impute it’s value into ˆy. 4 Prediction Models 4.1 Holt-Winter’s Holt-Winter’s algorithm [2] is used to forecast time series with trends and seasonal effects. In R, Holt-Winter’s algorithm is implemented in the stats package. 8 Friese et al. input : Univariate Time Series t input : Frequency fof t input : Constant kof similar windows to look for correlated values input : Method Mto calculate value to impute: mean, median, max or min output: Time Series ˆ t initialization; define a vector of window numbers to look for, using k win ← {−k, −k+ 1,...,−1,1,...,k−1, k} for i←initial point of data to end of data do if found ypoint is NA then look for correlated values to yand form vector tsby using win vector and freq. fas follows: ts← {y(tm−kf), y(tm−(k−1)f),...,y(tm−f), y(tm+f),...,y(tm+(k−1)f), y(tm+kf)} end switch method Mchosen, calculate do mean: ˆy←mean(ts) median: ˆy←median(ts) max: ˆy←max(ts) min: ˆy←min(ts) endsw end check if some NA values were imputed into time series ˆ t, and if there are, fix by using the left neighbor value. Algorithm 3: Imputation Algorithm. The method works with three exponential smoothing equations: Level: `t=α(yt−st−m) + (1 −α)(`t−1+bt−1) (4) Trend: bt=β∗(`t−lt−1+ (1 −β∗)bt−1(5) Seasonal: st=γ(yt−`t−1−bt−1) + (1 −γ)st−m(6) The forecasting equation is: Forecast: ˆyt+h|t=`t+hbt+st−m+h+ m(7) with h+ m=b(h−1) mod mc+ 1. α,β∗and γare so-called smoothing parameters. The level equation is a weighted average of the seasonally adjusted observation (yt−st−m) and the non-seasonal forecast (`t−1+bt−1) for time t. The trend equation shows that btis a weighted average of the estimated trend at time t based on `t−`t−1and bt−1, the previous estimate of the trend. The seasonal equation is a weighted average of the current seasonal index, (yt−`t−1−bt−1) and the seasonal index of the same season last term. Initial values for the level, trend, and seasonal indices are calculated using a simple decomposition and regression on the first two seasons of the given time series. The smoothing parameters are, if not given manually, fitted by the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method with the sum of squared errors of prediction (SSE). Univariate Frequency-based Imputation for Time Series Data 15 9 Summary and Outlook In this work, we introduced a scalable test problem generator for simulating missing data. We generated problem instances with small and large gaps, furthermore we introduced a new method called UniFIeD for handling missing data, which uses an automated frequency estimator using autocorrelation. We demonstrated that this method outperformes simple imputation methods. The results are independent from the underlying error measurements. Future plans involve the following steps: –extended experiments with additional data –discover the limits of the method –discover the maximum gap sizes –statistical validate the method –can it be used as a hybrid method –experiments on mutlivariate time-series We are also planning to provide UniFIeD as an R package on CRAN. 10 Acknowledgments This work has been kindly supported by the German Federal Ministry of Education and Research (BMBF) under the grants MCIOP (FKZ 17N0311) and CIMO (FKZ 17002X11). References 1. S. Buuren and K. Groothuis-Oudshoorn. Mice: Multivariate imputation by chained equations in r. Journal of Statistical Software, 45(3), 2011. 2. C. Holt. Forecasting seasonals and trends by exponentially weighted moving averages. International Journal of Forecasting, 20(1):5–10, 2004. 3. J. Honaker, G. King, and M. Blackwell. Amelia ii: A program for missing data. Journal of Statistical Software, 45(7):1–47, 12 2011. 4. N. Horton and K. Kleinman. Much ado about nothing. The American Statistician, 61(1):79–90, 2007. 5. R. Hyndman and Y. Khandakar. Automatic time series for forecasting: The forecast package for r. Technical report, Monash University, Department of Econometrics and Business Statistics, 2007. 6. R. J. Hyndman and Y. Khandakar. Automatic Time Series Forecasting. The forecast Package for R. Journal of Statistical Software, 27(3):1–22, 7 2008. 7. H. Junninen, H. Niska, K. Tuppurainen, J. Ruuskanen, and M. Kolehmainen. Methods for imputation of missing values in air quality data sets. Atmospheric Environment, 38(18):2895 – 2907, 2004. 8. R Development Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, 2011. ISBN 3-900051-07-0. 9. D. B. Rubin. Inference and missing data. Biometrika, 63(3):581–592, 1976. 10. A. Zeileis and G. Grothendieck. zoo: S3 infrastructure for regular and irregular time series. arXiv preprint math/0505527, 2005. Kontakt/Impressum Diese Ver¨offentlichungen erscheinen im Rahmen der Schriftenreihe ”CIplus”. Alle Ver¨offentlichungen dieser Reihe k¨onnen unter www.ciplus-research.de oder unter http://opus.bsz-bw.de/fhk/index.php?la=de abgerufen werden. K¨oln, Januar 2012 Herausgeber / Editorship Prof. Dr. Thomas Bartz-Beielstein, Prof. Dr. Wolfgang Konen, Prof. Dr. Horst Stenzel, Dr. Boris Naujoks Institute of Computer Science, Faculty of Computer Science and Engineering Science, Cologne University of Applied Sciences, Steinm¨ullerallee 1, 51643 Gummersbach url: www.ciplus-research.de Schriftleitung und Ansprechpartner/ Contact editor’s office Prof. Dr. Thomas Bartz-Beielstein, Institute of Computer Science, Faculty of Computer Science and Engineering Science, Cologne University of Applied Sciences, Steinm¨ullerallee 1, 51643 Gummersbach phone: +49 2261 8196 6391 url: http://www.gm.fh-koeln.de/~bartz/ eMail: [email protected] ISSN (online) 2194-2870