In this paper, we present a simple yet effective algorithm, called the Top-k Case Matching algorithm, for the imputation of missing values in streams of time series data that are similar to each other. The key idea of the algorithm is to look for the k situations in the historical data that are most similar to the current situation and to derive the missing value from the measured values at these k time points. To efficiently identify the top-k most similar historical situations, we adopt Fagin's Threshold Algorithm, yielding an algorithm with sub-linear runtime complexity with high probability, and linear complexity in the worst case (excluding the initial sorting of the data, which is done only once). We provide the results of a first experimental evaluation using real-world meteorological data. Our algorithm achieves a high accuracy and is more accurate and efficient than two more complex state of the art solutions.
Time series data is ubiquitous, e.g., in the financial stock market or in meteorology. In many applications time series data is incomplete, that is some values are missing for various reasons, e.g., sensor failures or transmission errors. However, many applications assume complete data, hence need to recover missing values before further data processing is possible.
In this paper, we focus on the imputation of missing values in long streams of meteorological time series data. As a case study, we use real-world meteorological data collected by the Südtiroler Beratungsring1 (SBR), which is an organization that provides professional and independent consultancy to the local wine and apple farmers, e.g., to determine the optimal harvesting time or to warn about potential threats, such as apple scab, fire blight, or frost. Es1http://www.beratungsring.org/ pecially frost is dangerous as it can destroy the harvest within a few minutes unless the farmers react immediately. The Südtiroler Beratungsring operates more than 120 weather stations spread all over South Tyrol, where each of them collects every five minutes up to 20 measurements including temperature, humidity etc. The weather stations frequently suffer outages due to sensor failures or errors in the transmission of the data. However, the continuous monitoring of the current weather condition is crucial to immediately warn about imminent threats such as frost and therefore the need arises to recover those missing values as soon as they are detected.
In this paper, we propose an accurate and efficient method to automatically recover missing values. The need for a continuous monitoring of the weather condition at the SBR has two important implications for our solution. Firstly, the proposed algorithm has to be efficient enough to complete the imputation before the next set of measurements arrive in a few minutes time. Secondly, the algorithm cannot use future measurements which would facilitate the imputation, since they are not yet available.
The key idea of our Top-k Case Matching algorithm is to seek for the k time points in the historical data when the measured values at a set of reference stations were most similar to the measured values at the current time point (i.e., the time point when a value is missing). The missing value is then derived from the values at the k past time points. While a naïve solution to identify the top-k most similar historical situations would have to scan the entire data set, we adopt Fagin’s Threshold Algorithm, which efficiently answers top-k queries by scanning, on average, only a small portion of the data. The runtime complexity of our solution is derived from the Threshold Algorithm and is sub-linear with high probability and linear in the worst case, when all data need to be scanned. We provide the results of a first experimental evaluation using real-world meteorological data from the SBR. The results are promising both in terms of efficiency and accuracy. Our algorithm achieves a high accuracy and is more accurate than two state of the art solutions.
The rest of the paper is organized as follows. In Section 2, we review the existing literature about imputation methods for missing values. In Section 3, we introduce the basic notation and a running example. In Section 4, we present our Top-k Case Matching algorithm for the imputation of missing values, followed by the results of an experimental evaluation in Section 5. Section 6 concludes the paper and outlines ideas for future work. 2.
Khayati et al. [
Khayati et al. [
The above techniques are designed to handle missing values in static time series. Therefore, they are not applicable in our scenario, as we have to continuously impute missing values as soon as they appear. A naïve approach to run the algorithms each time a missing value occurs is not feasible due to their relatively high runtime complexity.
There are numerous statistical approaches for the imputation of
missing values, including easy ones such as linear or spline
interpolation, all the way up to more complex models such as the ARIMA
model. The ARIMA model [
Let S = fs1; : : : ; sng be a set of time series. Each time series, s 2 S, has associated a set of reference time series Rs, Rs s . The value of a time series s 2 S at time t is denoted as S n f g s(t). A sliding window of a time series s is denoted as s([t1; t2]) and represents all values between t1 and t2.
EXAMPLE 1. Table 1 shows four temperature time series in a time window w = [1; 7], which in our application corresponds to seven timestamps in a range of 30 minutes. s is the base time series from the weather station in Schlanders, and Rs = fr1; r2; r3g is the associated set of reference time series containing the stations of Kortsch, Göflan, and Laas, respectively. The temperature value s(7) is missing. Figure 1 visualizes this example graphically.
The Top-k Case Matching algorithm we propose assumes that the time series data is aligned, which generally is not the case for our data. Each weather station collects roughly every 5 minutes new measurements and transmits them to a central server. Since the stations are not perfectly synchronized, the timestamps of the measurements typically differ, e.g., one station collects measurements at 09:02, 09:07, . . . , while another station collects them at 09:04, 09:09, . . . . Therefore, in a pre-processing step we align the time series data using linear interpolation, which yields measurement values every 5 minutes (e.g., 00:00, 00:05, 00:10, . . . ). If we observe a gap of more than 10 minutes in the measurements, we assume that the value is missing.
For the imputation of missing values we assign to each time series s a set Rs of reference time series, which are similar to s. The notion of similarity between two time series is tricky, though. Intuitively, we want time series to be similar when they have similar values and behave similarly, i.e., values increase and decrease roughly at the same time and by the same amount.
As a simple heuristic for time series similarity, we use the spatial proximity between the stations that record the respective time series. The underlying assumption is that, if the weather stations are nearby (say within a radius of 5 kilometers), the measured values should be similar, too. Based on this assumption, we manually compiled a list of 3–5 reference time series for each time series. This heuristic turned out to work well in most cases, though there are situations where the assumption simply does not hold. One reason for the generally good results is most likely that in our data set the over 100 weather stations cover a relatively small area, and hence the stations are very close to each other. 4.
Weather phenomena are often repeating, meaning that for example during a hot summer day in 2014 the temperature measured at the various weather stations are about the same as those measured during an equally hot summer day in 2011. We use this observation for the imputation of missing values. Let s be a time series where the current measurement at time , s( ), is missing. Our assumption on which we base the imputation is as follows: if we find historical situations in the reference time series Rs such that the past values are very close to the current values at time , then also the past measurements in s should be very similar to the missing value s( ). Based on this assumption, the algorithm searches for similar climatic situations in historical measurements, thereby leveraging the vast history of weather records collected by the SBR.
More formally, given a base time series s with reference time series Rs, we are looking for the k timestamps (i.e., historical situations), D = ft1; : : : ; tkg, ti < function , which minimize the error (t) = X jr( ) r2Rs r(t)j: That is, (t) (t0) for all t 2 D and t0 62 D [ f g. The error function (t) is the accumulated absolute difference between the current temperature r( ) and the temperature at time t, r(t), over all reference time series r 2 Rs. Once D is determined, the missing value is recovered using some aggregation function g (fs(t)j8t 2 Dg) over the measured values of the time series s at the timestamps in D. In our experiments we tested the average and the median as aggregation function (cf. Section 5).
EXAMPLE 2. We show the imputation of the missing value s(7) in Table 1 using as aggregation function g the average. For the imputation, we seek the k = 2 most similar historical situations. The two timestamps D = f4; 1g minimize (t) with (4) = j15:0° 15:0°j + j16:0° 15:9°j + j14:3° 14:2°j = 0:2° and (1) = 0:3°. The imputation is then simply the average of the base station measurements at time t = 4 and t = 1, i.e.,s(7) = avg(16:2°; 16:1°) = 21 (16:2° + 16:1°) = 16:15°.
A naïve implementation of this algorithm would have to scan the entire database of historical data to find the k timestamps that minimize (t). This is, however, not scalable for huge time series data, hence a more efficient technique is needed. 4.1
What we are actually trying to do is to answer a top-k query for
the k timestamps which minimize (t). There exist efficient
algorithms for top-k queries. For example, Fagin’s algorithm [
Our Top-k Case Matching algorithm is a variation of TA with slightly different settings. Fagin et al. assume objects with m attributes, a grade for each attribute and a monotone aggregation function f : Rm 7! R, which aggregates the m grades of an object into an overall grade. The monotonicity property is defined as follows.
DEFINITION 1. (Monotonicity) Let x1; : : : ; xm and x01; : : : ; x0m be the m grades for objects X and X0, respectively. The aggregation function f is monotone if f (x1; : : : ; xm) f (x01; : : : ; x0m) given that xi x0i for each 1 i m.
The TA finds the k objects that maximize the function f . To do so it requires two modes of accessing the data, one being sorted and the other random access. The sorted access is ensured by maintaining a sorted list Li for each attribute mi, ordered by the grade in descending order. TA keeps a bounded buffer of size k and scans each list Li in parallel until the buffer contains k objects and the lowest ranked object in the buffer has an aggregated grade that is greater than or equal to some threshold . The threshold is computed using the aggregation function f over the grades last seen under the sorted access for each list Li.
EXAMPLE 3. Table 2 shows four objects fA; B; C; Dg and their grade for the two attributes interestingness and popularity. Let us assume that k = 2 and the aggregation function f (x1; x2) = x1 + x2. Further, assume that the bounded buffer currently contains f(C; 18); (A; 16)g and the algorithm has read the data up to the boxes shown in gray. At this point the algorithm computes the threshold using the interestingness grade for object B and the popularity grade of object C, yielding = f (5; 9) = 5 + 9 = 14. Since the lowest ranked object in the buffer, object A, has an aggregated grade that is greater than , we can conclude that C and A are the top-2 objects. Note that the algorithm never read object D, yet it can conclude that D cannot be part of the top-k list.
interestingness
popularity
D
In order to use the Threshold Algorithm for the imputation of missing values in time series data, we have to adapt it. Instead of looking for the top-k objects that maximize the aggregation function f , we want to find the top-k timestamps that minimize the error function (t) over the reference time series Rs. Similar to TA, we need sorted access to the data. Therefore, for each time series r 2 Rs we define Lr to be the time series r ordered first by value and then by timestamp in ascending order. Table 3 shows the sorted data for the three reference time series of our running example (ignore the gray boxes and small subscript numbers for the moment).
The general idea of our modified TA algorithm is the following. The scan of each sorted lists starts at the current element, i.e., the element with the timestamp t = . Instead of scanning the lists Lri only in one direction as TA does, we scan each list sequentially in two directions. Hence, as an initialization step, the algorithm places two pointers, posr+ and posr , at the current value r( ) of time series r (the gray boxes in Table 3). During the execution of the algorithm, pointer posr+ is only incremented (i.e., moved down the list), whereas posr is only decremented (i.e., moved up the list). To maintain the k highest ranking timestamps, the algorithm uses a bounded buffer of size k. A new timestamp t0 is added only if the buffer is either not yet full or (t0) < (t), where t is the last (i.e., lowest ranking) timestamp in the buffer. In the latter case the timestamp t is removed from the buffer.
After this initialization, the algorithm iterates over the lists Lr in round robin fashion, i.e., once the last list is reached, the algorithm wraps around and continues again with the first list. In each iteration, exactly one list Lr is processed, and either pointer posr+ or posr is advanced, depending on which value the two pointers point to has a smaller absolute difference to the current value at time , r( ). This process grows a neighborhood around the element r( ) in each list. Whenever a pointer is advanced by one position, the timestamp t at the new position is processed. At this point, the algorithm needs random access to the values r(t) in each list to compute the error function (t). Time t is added to the bounded buffer using the semantics described above.
The algorithm terminates once the error at the lowest ranking timestamp, t, among the k timestamps in the buffer is less or equal to the threshold, i.e., (t) . The threshold is defined as = Pr2Rs jr( ) r(posr)j, where posr is either posr+ or posr , depending on which pointer was advanced last. That is, is the sum over all lists Lr of the absolute differences between r( ) and the value under posr+ or posr .
EXAMPLE 4. We illustrate the Top-k Case Matching algorithm for k = 2 and = 7. Table 4 shows the state of the algorithm in each iteration i. The first column shows an iteration counter i, the second the buffer with the k current best timestamps, and the last column the threshold . The buffer entries are tuples of the form (t; (t)). In iteration i = 1, the algorithm moves the pointer to t = 4 in list Lr1 and adds (t = 4; (4) = 0:2°) to the buffer. Since (4) = 0:2° > 0:0° = , the algorithm continues. The pointer in Lr2 is moved to t = 6, and (6; 0:4°) is added to the buffer. In iteration i = 4, timestamp 6 is replaced by timestamp 1. Finally, in iteration i = 6, the error at timestamp t = 1 is smaller or equal to , i.e., (1) = 0:3° 6 = 0:3°. The algorithm terminates and returns the two timestamps D = f4; 1g.
Iteration i 1 2 3 4 5 6
Buffer (4; 0:2°) (4; 0:2°); (6; 0:4°) (4; 0:2°); (6; 0:4°) (4; 0:2°); (1; 0:3°) (4; 0:2°); (1; 0:3°) (4; 0:2°); (1; 0:3°)
Threshold i 0:0° 0:0° 0:1° 0:1° 0:2° 0:3°
Algorithm 1 shows the pseudo code of the Top-k Case Matching algorithm. The algorithm has three input parameters: a set of time series Rs, the current timestamp , and the parameter k. It returns the top-k most similar timestamps to the current timestamp . In line 2 the algorithm initializes the bounded buffer of size k, and in line 4 the pointers posr+ and posr are initialized for each reference time series r 2 Rs. In each iteration of the loop in line 7, the algorithm advances either posr+ or posr (by calling Algorithm 2) and reads a new timestamp t. The timestamp t is added to the bounded buffer using the semantics described before. In line 15, the algorithm computes the threshold . If the buffer contains k timestamps and we have (t) , the top-k most similar timestamps were found and the algorithm terminates.
Algorithm 2 is responsible for moving the pointers posr+ and posr for each list Lr. The algorithm uses three utility functions. The first is next(), which takes a pointer as input and returns the next position by either incrementing or decrementing, depending
on the direction of the pointer. If next() reaches the end of a list, it returns NIL. The utility functions timestamp() and value() return the timestamp and value of a list Lr at a given position, respectively. There are four cases, which the algorithm has to distinguish: 1. None of the two pointers reached the beginning or end of the list. In this case, the algorithm checks which pointer to advance (line 5). The pointer that is closer to r( ) after advancing is moved by one position. In case of a tie, we arbitrarily decided to advance posr+. 2. Only posr reached the beginning of the list: the algorithm increments posr+ (line 11). 3. Only posr+ reached the end of the list: the algorithm decrements posr (line 13). 4. The two pointers reached the beginning respective end of the list: no pointer is moved.
In the first three cases, the algorithm returns the timestamp that was discovered after advancing the pointer. In the last case, NIL is returned.
At the moment we use an in-memory implementation of the algorithm, which loads the whole data set into main memory. More specifically, we keep two copies of the data in memory: the data sorted by timestamp for fast random access and the data sorted by value and timestamp for fast sorted access.
Note that we did not normalize the raw data using some standard technique like the z-score normalization, as we cannot compute that efficiently for streams of data without increasing the complexity of our algorithm. 4.4
The correctness of the Top-k Case Matching algorithm follows directly from the correctness of the Threshold Algorithm. What remains to be shown, however, is that the aggregation function (t) is monotone.
Result: Next timestamp to look at or NIL
THEOREM 4.1. The aggregation function (t) is a monotonically increasing function.
PROOF. Let t1 and t2 be two timestamps such that jr( ) r(t1)j jr( ) r(t2)j for each r 2 Rs. Then it trivially follows that (t1) (t2) as the aggregation function is the sum of jr( ) r(t1)j over each r 2 Rs and, by definition, each component of (t1) is less than or equal to the corresponding component in (t2). 4.5
The space and runtime bounds of the algorithm follow directly from the probabilistic guarantees of TA, which has sub-linear cost with high probability and linear cost in the worst case. Note that sorting the raw data to build the lists Lr is a one-time preprocessing step with complexity O(n log n). After that the system can insert new measurements efficiently into the sorted lists with logarithmic cost.
In this section, we present preliminary results of an experimental
evaluation of the proposed Top-k Case Matching algorithm. First,
we study the impact of parameter k on the Top-k Case Matching
and a baseline algorithm. The baseline algorithm, referred to as
“Simple Average”, imputes the missing value s( ) with the average
of the values in the reference time series at time , i.e., s( ) =
jR1sj Pr2Rs r( ). Second, we compare our solution with two state
of the art competitors, REBOM [
In this experiment, we study the impact of parameter k on the accuracy and the runtime of our algorithm. We picked five base stations distributed all over South Tyrol, each having two to five reference stations. We simulated a failure of the base station during a time interval, w, of 8 days in the month of April 2013. This amounts to a total of 11452 missing values. We then used the Top-k Case Matching (using both the average and median as aggregation function g) and Simple Average algorithms to impute the missing values. As a measure of accuracy we use the average absolute difference between the real value s( ) and the imputed value s ( ), i.e., = jw1 j P 2w js( ) s ( )j
Figure 2 shows how the accuracy of the algorithms changes with varying k. Interestingly and somewhat unexpectedly, decreases as k increases. This is somehow contrary to what we expected, since with an increasing k also the error function (t) grows, and therefore less similar historical situations are used for the imputation. However, after a careful analysis of the results it turned out that for low values of k the algorithm is more sensitive to outliers, and due to the often low quality of the raw data the imputation is flawed. 50
Parameter k
Table 5 shows an example of flawed raw data. The first row is the current situation, and we assume that the value in the gray box is missing and need to be recovered. The search for the k = 3 most similar situations using our algorithm yields the three rows at the bottom. Notice that one base station value is 39:9° around midnight of a day in August, which is obviously a very unlikely thing to happen. By increasing k, the impact of such outliers is reduced and hence decreases. Furthermore, using the median as aggregation function reduces the impact of outliers and therefore yields better results than the average.
Timestamp
Figure 3 shows the runtime, which for the Top-k Case Matching algorithm linearly increases with k. Notice that, although the imputation of missing values for 8 days takes several minutes, the algorithm is fast enough to continuously impute missing values in our application at the SBR. The experiment essentially corresponds to a scenario, where in 11452 base stations an error occurs at the same time. With 120 weather stations operated by the SBR, the number of missing values at each time is only a tiny fraction of the missing values that we simulated in this experiment. 5.2
In this experiment, we compare the Top-k Case Matching
algorithm with two state-of-the-art algorithms, REBOM [
Parameter k
The box plot in Figure 4 shows how the imputation error js( ) s ( )j is distributed for each of the four algorithms. The left and right line of the box are the first and third quartile, respectively. The line inside the box denotes the median and the left and right whiskers are the 2.5% and 97.5% percentile, which means that the plot incorporates 95% of the values and omits statistical outliers. The experiment clearly shows that the Top-k Case Matching algorithm is able to impute the missing values more accurately than CD and REBOM. Although not visualized, also the maximum observed error for our algorithm is with 2:29° (Average) and 2:21° (Median) considerably lower than 3:71° for CD and 3:6° for REBOM.
0
CD REBOM 0
0:5 1 1:5 Absolute Difference in °C 2
In terms of runtime, the Top-k Case Matching algorithm needed 16 seconds for the imputation of the 2017 missing measurements, whereas CD and REBOM needed roughly 10 minutes each. Note, however, that this large difference in run time is also due to the fact that CD and REBOM need to compute the Pearson correlation coefficient which is a time intensive operation.
In this paper, we presented a simple yet efficient and accurate algorithm, termed Top-k Case Matching, for the imputation of missing values in time series data, where the time series are similar to each other. The basic idea of the algorithm is to look for the k situations in the historical data that are most similar to the current situation and to derive the missing values from the data at these time points. Our Top-k Case Matching algorithm is based on Fagin’s Threshold Algorithm. We presented the results of a first experimental evaluation. The Top-k Case Matching algorithm achieves a high accuracy and outperforms two state of the art solutions both in terms of accuracy and runtime.
As next steps we will continue with the evaluation of the
algorithm, taking into account also model based techniques such as
DynaMMo [
There are several possibilities to further improve the algorithm. First, we would like to explore whether the algorithm can dynamically determine an optimal value for the parameter k, which is currently given by the user. Second, we would like to make the algorithm more robust against outliers. For example, the algorithm could consider only historical situations that occur roughly at the same time of the day. Moreover, we can bend the definition of “current situation” to not only consider the current timestamp, but rather a small window of consecutive timestamps. This should make the ranking more robust against anomalies in the raw data and weather phenomena such as the foehn. Third, right now the similarity between time series is based solely on temperature data. We would like to include the other time series data collected by the weather stations, such as humidity, precipitation, wind, etc. Finally, the algorithm should be able to automatically choose the currently handpicked reference time series based on some similarity measures, such as the Pearson correlation coefficient.
The work has been done as part of the DASA project, which is funded by the Foundation of the Free University of Bozen-Bolzano. We wish to thank our partners at the Südtiroler Beratungsring and the Research Centre for Agriculture and Forestry Laimburg for the good collaboration and helpful domain insights they provided, in particular Armin Hofer, Martin Thalheimer, and Robert Wiedmer. 8.