Repository: sryza/spark-timeseries (language: python)

Issue:
Tick frequency?
I notice that increasingly we end up having to implement algorithms for a sequence and then the same algorithm for a time step. In theory a time series is something that is based on time, not on indexing (like a sequence is), but we seem to be using both more or less interchangeably.

Maybe if we added a Frequency, say TickFrequency, that represents indexing rather than a notion of time, we could unify each version of each algorithm (lagging, differencing, etc.) into one, thus simplifying the API?

Lagging irregular time series
I am now faced with the problem of regressing on lags of a time series that is irregular. So next I will tackle the problem of lagging an irregular time series. I'd be curious to see if anyone has ideas or opinions on how to handle this problem.

A high-level overview of my proposal at the moment is as follows.

1 - The caller of the function specifies the base frequency to use. We can later add a function that automatically discovers the base frequency by looking at the minimal time gap that exists in the time series. But it's good in any case to have a manual override, so I will start there.

2 - The caller will provide an interpolation function that tells us what to do when there is no value at time (t - lag). This function will have the form (prev, next) => Double, where prev and next will represent the values surrounding time (t - lag) on the lower and upper boundaries. 

3 - The resulting lags will look something like this if we specify lag periods of 0.1 and 0.2:

Time, Value
0.1, 100
0.2, 200
0.24, 123
0.25, 300
0.4, 400
0.5, 500

becomes:

Time, Lag0.1(Value), Lag0.2(Value)
0.3, 200, 100
0.4, f(300, 400), 200
0.5, 400, f(300, 400)

where of course f is the user-provided interpolation function. 

Notice that the time instant 0.3 gets created even though it didn't exist in the previous dataset. Also, value (0.24, 123) gets completely ignored: there is loss of information because the specified base frequency is lower than the "true" base frequency of the dataset.



Discussion / hints:
Would you mind clarifying what you mean by time step and sequence?  Are there cases where it's not possible to resample a time series, so that these line up?

A sequence is an ordered list of values. The differences() currently implemented in TimeSeries, for example, is a difference on sequences. It doesn't care about time.

A time series is when each observation is associated to a timestamp. In this case, we may want to difference every hour, for example. If the time series is irregular (and yes, if the time series is uniform, the time notion can be dropped and it can be treated as a sequence), then lagging or differencing by an hour can mean a different index delta for each observation.

It's always (I think) possible to resample and interpolate an irregular times series to turn it into a uniform one -- but to me the right way to handle this is to treat the "tick" frequency (i.e. sequence indexes) as just another time frequency. 

I used to think that irregular time series are fundamentally useless and bound to be transformed into uniform eventually in order to do something with them -- but I'm less and less convinced of that. In fact, I now tend to keep my irregular time series irregular (thus I wouldn't want a library that forces us to transform it to uniform in order to perform simple lag/difference operations). Not to mention possible performance and memory consumption inefficiencies from doing this kind of two-step approach.

But I'm not sure what specific point you take issue with: you're saying we should only have sequence-based operations and have the user transform their time series to a uniform one all the time (and then calculate how many index positions correspond to the time period he's interested in)? 

Apologies, but I think I'm having trouble understanding how this would concretely play out in the API.  What would be the semantics of the TickFrequency implementation of Frequency's advance and difference methods?  Would the signature of TimeSeries.differences change?

Yes, the signature of differences would change, and would now take in a frequency in addition to the "lag" parameter. If the provided frequency is the tick frequency, then the semantic of the difference function would be the same as the current one (although the implementation would be more generic to support other frequencies).

The semantic of the TickFrequency's advance and difference methods would be trivial: advance adds one to the index, and difference subtracts the two provided indices.

The behaviour of differences for non-tick frequencies would probably be along the lines of taking the difference between each time t and the last value at or before time t - lag. The result would be an irregular index if the input was irregular (i.e. no interpolation or "uniformisation").

The advantage from a standpoint of code maintenance, redundancy elimination and API simplicity is that we would now have only 1 function to do differences (and 1 function to do lags, etc.)

> The semantic of the TickFrequency's advance and difference methods would be trivial: advance adds one to the index, and difference subtracts the two provided indices.

`Frequency.advance` accepts and returns a `ZonedDateTime`, so what would be the meaning of "adding one"?
Similarly, `Frequency.difference` accepts two `ZonedDateTime`s, so what would be the meaning of "subtracting" them?

Well I didn't express myself correctly. What I mean is that Frequency.advance(dt, n) will add n "ticks" (i.e. index positions) to dt. So essentially the implementation would be to look up the index at dt, add n, and reverse lookup the corresponding datetime.

Similarly for Frequency.difference, it looks up the two indices of the two provided datetimes, and returns the number of index positions between. It's equivalent to asking "How many trades between 8:05 and 8:35" for a financial context, for example.

What this means for an algorithm such as differences(n), is that we could have something along the following lines (to "pseudo-code" it, and ignoring some implementation details):

differences(n: Int, freq: Frequency):
{
  startLoc = locAtOrAfterDateTime(freq.advance(index.start(), n))

  difference vectors = map(
    for each entry from startLoc to the end of the time series
    {
      previous index = locAtOrBeforeDateTime(freq.advance(datetimeAtLoc(currentLoc), -n))
      difference at current index = current value - previous value
    }
  )

  new index = index.islice(startLoc, index.size)
  new TimeSeries(difference vectors, new index)
}

and voila, we have a single algorithm and function that works for both index-based and time-based "arithmetic". I'm sure we can come up with something equivalent for irregular lags. 

This behavior preserves the irregular index, rather than "uniformizing" it. 

Also, it's a bit sub-optimal in terms of performance for the tick case (as is often the case for more generic solutions), because it adds extra array lookups per entry. If we care, we can add an if statement that checks if it's a TickFrequency and avoids extra lookups (we still get 1 coherent API function for both cases, although more cyclomatically complex -- it's a matter of which tradeoff we want).

Does that make more sense?

I see.  So the TickFrequency would wrap a DateTimeIndex,.  It seems to me that this is basically circumventing the meaning of a Frequency as something that helps map between date-times and indices.  It also seems weird that the DateTime would have a reference to the Frequency, and the Frequency would have a reference to the DateTimeIndex.  Lastly, if we have a DateTimeIndex with observations every day, and the user wants to calculate differences for two day intervals - how would the user specify whether they want the odd days or the even days?

So my feeling is that the API would be more elegant to not offer a differences function that takes a Frequency at all, and instead require users to resample to a new DateTimeIndex with the desired Frequency. 

I completely forgot the part that Frequencies don't have access to the datetime array... You're right, that makes things a little too weird.

That being said, in my use case I need to be able to obtain differences without resampling (i.e. I want to preserve the irregular index). Maybe I'm misunderstanding what you mean by resampling.

So, to be more concrete, the behavior I need for irregular lags is as follows:

t, value
0.1, 3.5
0.3, 3.6
1.0, 3.7
1.1, 3.8
1.5, 3.9
2.1, 4.0
2.3, 4.1

lagged by 1 unit should give:

t, value
1.1, 3.5
1.5, 3.6
2.1, 3.8
2.3, 3.8

I'm not sure I understand how resampling solves this problem

I see.

So the way I think that resampling would address this would be by resampling to .1 granularity, with "fill with next".  So after resampling you'd end up with:

0.1, 3.5
0.2, 3.6
0.3, 3.6
0.4, 3.7
0.5, 3.7
0.6, 3.7
0.7, 3.7
0.8, 3.7
0.9, 3.7
1.0, 3.7
1.1, 3.8
1.2, 3.9 
1.3, 3.9
1.4, 3.9
1.5, 3.9
etc.

And then one could call `lag(10)`.

Would that work? Is your concern that the greatest common factor of the numbers could be too small, so it would blow up the time series?  If that's the issue, I wonder if we could introduce a new uniform DateTimeIndex implementation that represents the above time series in a sparse way.  I.e. it would wrap an IrregularDateTimeIndex, but provide values at uniform intervals.  If asked for the value at time 0.2, if it didn't contain a value for 0.2, it would find the next date-time that did have a value and return that.

My biggest concern is with downsampling -- due to the information loss it causes (for my use case).

In the case you propose, i.e. upsampling, I'm definitely a bit concerned with blowing up the time series. Initially that's what I was doing, but then I hit cases where the time series' precision was microsecond, and realized that if I upsampled it entirely to the microsecond it would produce an enormous time series to analyze. Bottom line, it potentially because a huge performance problem.

And, by the way, there is extra trickiness to zero-order hold upsampling like that: it only makes sense if the time series is I(1). If it's an I(0) time series, you must fill it in 0, not with the last value. For example: stock prices versus returns. You upsample prices by keeping the last value, but if you upsample returns by keeping the last value, it takes a whole different meaning (i.e. that the stock price keeps rising).

I'm really convinced that we need an algorithm that works in irregular indices and preserves them, rather than manipulating the data itself (which is always perilous). I have local versions that are simply extra functions, timeDifferences() and irregularLag(), but I didn't find that very elegant (due to 2 different versions of each functionality). The TickFrequency thing was my failed attempt at providing a more coherent and uniform API.

Regarding what you propose (the wrapped IrregularDateTimeIndex), I'm not sure yet. I need to think about it.

> In the case you propose, i.e. upsampling, I'm definitely a bit concerned with blowing up the time series...

Yeah, I agree this is definitely a concern that I hadn't considered earlier.

> I'm really convinced that we need an algorithm that works in irregular indices and preserves them, rather than manipulating the data itself (which is always perilous).

One of the main reasons I'm arguing for the resampling approach is that I think the alternative requires each algorithm to accept a ballooning amount of information.  For example, consider your point about how we should impute missing values for an I(0) time series differently than for an I(1) time series.  It's not difficult to imagine a similar example to the one you provided a couple days ago where one would want to impute 0s instead of the next present value in the irregular series.  If we wanted to support both imputation strategies, we'd need to add an extra argument to the lag function, as well as to all the other functions that you were suggesting should accept a TickFrequency.  Adding the imputation strategy instead to a resample function makes things more composable .

The signature of the lags function for Irregular datetime indexes should probably look like:

lags[U](maxLag: Int, lagPeriod: Period, includeOriginals: Boolean, laggedKey: %28K, Int%29 => U = laggedPairKey[K]_, interpolation: %28Double, Double%29 => Double): TimeSeries[U]

I will later create an equivalent for the lagsPerColumn as well.

Note: the date/time arithmetic to fetch the lags will have to be DST-aware (thus use the nscala-time concepts that are DST-aware).

Correction for the above post (the function signature didn't show up correctly):

lags[U: ClassTag](
      maxLag: Int,
      lagFrequency: Frequency,
      includeOriginals: Boolean,
      laggedKey: (K, Int) => U = laggedPairKey[K]_,
      interpolation: (Double, Double) => Double)
    : TimeSeries[U]

@sryza let me know if you have any thoughts on this

Although now I'm still facing that annoying

"Error:(24, 7) in class TimeSeries, multiple overloaded alternatives of method lags define default arguments.
class TimeSeries[K](val index: DateTimeIndex, val data: DenseMatrix[Double],
      ^"

compilation error... looks like I may have to give this function a unique name, such as lagIrregular

Signature for the "perColumn" version:

lagsPerColumnIrregular[U: ClassTag](
lagsPerCol: Map[K, (Boolean, Int)],
lagFrequency: Frequency,
includeOriginals: Boolean,
laggedKey: (K, Int) => U = laggedPairKey[K]_,
interpolation: (Double, Double) => Double)
: TimeSeries[U]

I'm realizing the interpolation function signature I proposed above is too narrow to accomodate most useful cases. All it can accomplish is zero-order hold and averaging of the 2 values...

A more useful linear interpolation would need to know about the datetimes corresponding to the 2 boundary values. But even so, most non-trivial interpolation methods require much more than the 2 neighbouring points.

In order to resolve this annoying signature issue, I suggest removing the default value from the original `lags` method

Original `lags` method

``` scala
def lags[U: ClassTag](maxLag: Int, includeOriginals: Boolean,
      laggedKey: (K, Int) => U = laggedPairKey[K]_)
    : TimeSeries[U] = {
```

`lags` method after removing the default value

``` scala
def lags[U: ClassTag](maxLag: Int, includeOriginals: Boolean,
 laggedKey: (K, Int) => U): TimeSeries[U] = {
```

and provide another `lags` overload without the `laggedKey` parameter that calls the original `lags` using `laggedPairKey`

``` scala
def lags[U: ClassTag](maxLag: Int, includeOriginals: Boolean): TimeSeries[U] = {
  lags[U](maxLag, includeOriginals, laggedPairKey[K])
}
```

How do you choose the first entry in the lagged index (e.g. 0.3 in your example)? I think that value affects the lagged time series considerably.

I can see that the lagged time series would have a uniform index. So, the irregular time series is implicitly being converted into a uniform time series. Why wouldn't we provide a `toUniform` utility that converts irregular indexes to uniform? and maybe irregular lags to use this utility internally?

I agree with @ahmed-mahran regarding uniform vs. irregular.  `lag` for an `IrregularDateTimeIndex` should just use the previous element in the sequence of datetimes composing the index, i.e. it shouldn't implicitly overlay a uniform index.  But there should be an easy way to get from an irregular index to a uniform index if that's the user's desire.

Regarding the uniform vs irregular -- fair enough. I suppose I'll call my function toUniform(), as this is the current use case I have. It's true that in a sense what I'm describing is more than just lagging, it's also converting to a uniform time series (explicitly, in fact, because when I instantiate the TimeSeries to be returned, I pass in a new UniformDatetimeIndex instance).

As for lagging irregular time series in an irregular way, I'm not clear what is the "meaning" of that. What does it really "mean" to lag by 0.4 seconds just because the previous tick was 0.4 seconds ago, and then by 1 second because the second previous tick was 1 second ago. I have a hard time imagining a context in which that makes sense, to be honest. Do you guys have any examples?

To answer @ahmed-mahran regarding choosing the first value (0.3), what I do is take the minimum (most anterior) time instant that exists in the series, and I add to it the maximum lag period that is specified by the user (so 0.1 + 0.2 in my example). Does that make sense? Do you have alternative approaches in mind that could work better?

BTW there's a mistake in my perColumn signature function above: the includeOriginals parameter should be removed.

Another important point I just realized: the interpolation function should be specific to each individual variable within the multivariate time series. 

For example, if my 2 individual time series are for the price of a stock and the volume traded for that stock, then the interpolation for the price might be zero-order hold (i.e. last traded price), but the interpolation for the volume would probably be a summation of the values since the last regularized tick (this is especially relevant if we are downsampling the time series)...

Regarding lagging irregular ts in irregular way, I think the nearest thing to @sryza's idea is to add the lag value to each instant with which to expand the time index. However, that would lead to a collection of sparse time series with a lot of NA's. Applying this on your example,

if we specify lag periods of 0.1 and 0.2:

| Time | Value |
| --- | --- |
| 0.1 | 100 |
| 0.2 | 200 |
| 0.24 | 123 |
| 0.25 | 300 |
| 0.4 | 400 |
| 0.5 | 500 |

becomes:

| Time | Lag0.1(Value) | Lag0.2(Value) |
| --- | --- | --- |
| 0.2 | 100 | NA |
| 0.3 | 200 | 100 |
| 0.34 | 123 | NA |
| 0.35 | 300 | NA |
| 0.4 | NA | 200 |
| 0.44 | NA | 123 |
| 0.45 | NA | 300 |
| 0.5 | 400 | NA |
| 0.6 | 500 | 400 |
| 0.7 | NA | 500 |

After that the user is free to carry out NA imputation, maybe time series compression. Then toUniform if desired.

An alternative solution is call toUniform then lag.

What do you think guys; I think we should open another separate issue for converting irregular to uniform?


Fix patch:
diff --git a/src/main/scala/com/cloudera/sparkts/TimeSeries.scala b/src/main/scala/com/cloudera/sparkts/TimeSeries.scala
index 5d42cf51..6b6cb9a6 100644
--- a/src/main/scala/com/cloudera/sparkts/TimeSeries.scala
+++ b/src/main/scala/com/cloudera/sparkts/TimeSeries.scala
@@ -18,7 +18,7 @@ package com.cloudera.sparkts
 import java.time._
 
 import breeze.linalg.{Axis, diff, DenseMatrix => BDM, DenseVector => BDV}
-import org.apache.spark.mllib.linalg.{DenseMatrix, Vector}
+import org.apache.spark.mllib.linalg.{DenseMatrix, Vector, Vectors}
 import MatrixUtil._
 import TimeSeries._
 
@@ -168,6 +168,72 @@ class TimeSeries[K](val index: DateTimeIndex, val data: DenseMatrix,
     new TimeSeries[K](index, mat, keys :+ key)
   }
 
+  /**
+    * Returns a TimeSeries where each individual time series is differenced with the given time
+    * frequency.
+    *
+    * Example:
+    * t    value
+    * 1pm  3.5
+    * 2pm  3.6
+    * 10pm 4.6
+    * 11pm 5.9
+    * 12pm 6.6
+    * And we difference by a frequency of 10 hours, we get:
+    *
+    * t    d(value)
+    * 11pm 2.4
+    * 12pm 3.0
+    *
+    * The process is as follows: for each timestamp t, find the value at (t - frequency).
+    * If there is none, go back to the previous available value. If there is none, discard t
+    * because we cannot difference it. This results in discarding all values earlier than time
+    * after the first available value in the time series, because they cannot be differenced due
+    * to missing information.
+    *
+    * It handles the NaN values in the following way:
+    *
+    * If we are differencing [x(t) - x(t-k)], for example, where k is the lag period, then
+    * if x(t) is NaN, the result is NaN. If x(t-k) is NaN but x(t) is not, however, it looks for
+    * the most recent non-NaN value before t-k and uses that to calculate the difference.
+    */
+  def differencesByFrequency(baseFrequency: Frequency): TimeSeries[K] = {
+    // first calculate the new index
+    val startDateTime = baseFrequency.advance(index.first, 1)
+    var startIndex = index.locAtOrAfterDateTime(startDateTime)
+    if (startIndex == 0) {
+      startIndex = 1
+    }
+    val newIndex = index.islice(startIndex, index.size)
+
+    // now calculate the differences of the time series
+    mapSeries(vec => {
+      var output = Array[Double]()
+      for (i <- startIndex until index.size) {
+        if (vec(i).isNaN()) {
+          output = output :+ Double.NaN
+        } else {
+          val prevInstant = baseFrequency.advance(index.dateTimeAtLoc(i), -1)
+          var prevLoc = index.locAtOrBeforeDateTime(prevInstant)
+
+          while (vec(prevLoc).isNaN() && prevLoc > 0) {
+            prevLoc = prevLoc - 1
+          }
+          if (prevLoc == -1) prevLoc = 0
+
+          if (vec(prevLoc).isNaN()) {
+            output = output :+ Double.NaN // (any value - NaN) = NaN
+          } else {
+            val valueDiff = vec(i) - vec(prevLoc)
+            output = output :+ valueDiff
+          }
+        }
+      }
+
+      Vectors.dense(output)
+    }, newIndex)
+  }
+
   /**
    * Returns a TimeSeries where each time series is differenced with the given order. The new
    * TimeSeries will be missing the first n date-times.


Test patch:
diff --git a/src/test/scala/com/cloudera/sparkts/TimeSeriesSuite.scala b/src/test/scala/com/cloudera/sparkts/TimeSeriesSuite.scala
index ca44d91b..91b9f7a6 100644
--- a/src/test/scala/com/cloudera/sparkts/TimeSeriesSuite.scala
+++ b/src/test/scala/com/cloudera/sparkts/TimeSeriesSuite.scala
@@ -106,4 +106,71 @@ class TimeSeriesSuite extends FunSuite with ShouldMatchers {
     filteredTS.toInstants()(0)._1.isAfter(dt.plusDays(1)) should be (true)
     filteredTS.toInstants()(1)._1.isAfter(dt.plusDays(1)) should be (true)
   }
+
+  test("differencesByFrequency example 1") {
+    val dt = ZonedDateTime.of(2015, 4, 8, 0, 0, 0, 0, ZoneId.of("Z"))
+    val samples = Array(
+      ((dt, Array(1.0, 2.0, 3.0))),
+      ((dt.plusMinutes(1), Array(4.0, 5.0, 6.0))),
+      ((dt.plusMinutes(2), Array(5.5, 6.5, 7.0))),
+      ((dt.plusMinutes(4), Array(7.0, 8.0, 9.0))),
+      ((dt.plusMinutes(7), Array(7.0, 8.0, 9.0))),
+      ((dt.plusMinutes(7).plusSeconds(5), Array(7.1, 8.2, 9.3))),
+      ((dt.plusMinutes(7).plusSeconds(15), Array(7.4, 8.4, 9.4))),
+      ((dt.plusMinutes(8), Array(6.0, 7.0, 8.0))),
+      ((dt.plusMinutes(12), Array(7.0, 8.0, 9.0))),
+      ((dt.plusMinutes(14), Array(10.5, 11.6, 12.7)))
+    )
+
+    val labels = Array("a", "b", "c")
+    val ts = timeSeriesFromIrregularSamples(samples, labels)
+
+    val differencedTS = ts.differencesByFrequency(new MinuteFrequency(1))
+
+    differencedTS.toInstants()(0)._1.getMinute() should be (1)
+    differencedTS.toInstants()(1)._1.getMinute() should be (2)
+    differencedTS.toInstants()(2)._1.getMinute() should be (4)
+    differencedTS.toInstants()(3)._1.getMinute() should be (7)
+    differencedTS.toInstants()(4)._1.getMinute() should be (7)
+    differencedTS.toInstants()(5)._1.getMinute() should be (7)
+    differencedTS.toInstants()(6)._1.getMinute() should be (8)
+    differencedTS.toInstants()(7)._1.getMinute() should be (12)
+    differencedTS.toInstants()(8)._1.getMinute() should be (14)
+
+    differencedTS.toInstants()(0)._2.toArray should be (Array(3.0, 3.0, 3.0))
+    differencedTS.toInstants()(1)._2.toArray should be (Array(1.5, 1.5, 1.0))
+    differencedTS.toInstants()(2)._2.toArray should be (Array(1.5, 1.5, 2.0))
+    differencedTS.toInstants()(3)._2.toArray should be (Array(0.0, 0.0, 0.0))
+    arrayShouldBe(differencedTS.toInstants()(4)._2.toArray, Array(0.1, 0.2, 0.3), 0.01)
+    arrayShouldBe(differencedTS.toInstants()(5)._2.toArray, Array(0.4, 0.4, 0.4), 0.01)
+    arrayShouldBe(differencedTS.toInstants()(6)._2.toArray, Array(-1.0, -1.0, -1.0), 0.01)
+    arrayShouldBe(differencedTS.toInstants()(7)._2.toArray, Array(1.0, 1.0, 1.0), 0.01)
+    arrayShouldBe(differencedTS.toInstants()(8)._2.toArray, Array(3.5, 3.6, 3.7), 0.01)
+  }
+
+  test("differencesByFrequency example 2") {
+    val dt = ZonedDateTime.of(2015, 4, 8, 1, 0, 0, 0, ZoneId.of("Z"))
+    val samples = Array(
+      ((dt, Array(3.5))),
+      ((dt.plusHours(1), Array(3.6))),
+      ((dt.plusHours(9), Array(4.6))),
+      ((dt.plusHours(10), Array(5.9)))
+    )
+
+    val labels = Array("a")
+    val ts = timeSeriesFromIrregularSamples(samples, labels)
+
+    // This should only return 1 timestamp: 11pm, with value: 2.4
+    val differencedTS = ts.differencesByFrequency(new HourFrequency(10))
+
+    differencedTS.head._2.size shouldBe (1)
+    differencedTS.head._2(0) shouldBe (2.4 +- 0.01)
+  }
+
+  def arrayShouldBe(left: Array[Double], right: Array[Double], tol: Double) = {
+    for (i <- 0 until left.length)
+    {
+      left(i) should be (right(i) +- tol)
+    }
+  }
 }