I'm trying to adapt an algorithm to calculate covariance to work over a rolling window on the data. Wikipedia has an algorithm for online covariance:

def online_covariance(data1, data2):
    mean1 = mean2 = 0
    M12 = 0
    for x, y in zip(data1, data2):
        n += 1
        delta1 = (x - mean1) / n
        mean1 += delta1
        delta2 = (y - mean2) / n
        mean2 += delta2
        M12 += (n - 1) * delta1 * delta2 - M12 / n
    return n / (n - 1) * M12

But I need this to work over an arbitrarily sized rolling window. I've already adapted a different algorithm on that Wikipedia page for variance to work on a rolling window, but I'm getting stuck doing it for covariance.

The below pseudo-code is what I have so far. Assume the window function returns an enumerable over a structure that contains the elements entering and leaving the window (Entered and Exited, respectively).

def online_covariance(data1, data2, windowSize):
    mean1 = mean2 = 0
    M12 = 0
    i = 0
    covars = []
    for x, y in zip(window(data1, windowSize), window(data2, windowSize)):

        if(n < windowSize)
            n += 1
        delta1 = (i < windowSize ? x.Entered - mean1 : x.Entered - x.Exited) / n
        mean1 += delta1
        delta2 = (i < windowSize ? y.Entered - mean2 : y.Entered - y.Exited) / n
        mean2 += delta2

        // Up to this point I can verify the algorithm is correct- 
        // it keeps the correct means for each data set
        // This next line is where I'm stuck- 
        // it needs to be modified to remove the value that 
        // just left the window
        M12 += (n - 1) * delta1 * delta2 - M12 / n

        covars.Add(n / (n - 1) * M12)
    return covars
  • Have you tried just passing the data for your rolling window into the original function? – Robert Harvey Dec 6 '16 at 17:08
  • @RobertHarvey I need to calculate the results of consecutive rolling windows over a large date range. Calling the function over and over again would be very slow- thus wanting to do this with an online algorithm. – MgSam Dec 6 '16 at 17:54
  • Well, couldn't you pass your data range by reference, and the start and end points, to your function? You could then loop over that, outside of the function. Calling the function repeatedly is only going to be slower if you're copying the data in your window and passing the copy to the function. – Robert Harvey Dec 6 '16 at 17:58
  • I'm saying that you're having difficulty because your function is doing too much. Make a function that computes one window given the start and end locations within your large data set, and then make another function that loops over those windows. – Robert Harvey Dec 6 '16 at 17:59
  • 1
    But don't you still have to loop over the subset anyway? Or is it a matter of keeping a running total? I haven't looked at the problem in detail, but couldn't you simply subtract the value that just left the window, and add the value that just entered the window? – Robert Harvey Dec 6 '16 at 18:06

To begin with

Just think about how you'd compute Var(X) online (where X is an unending stream of numbers).

What you'd do here is keep track of two quantities:

  1. Sum of all X values (i.e. sum(x1, x2, x3, ....))
  2. Sum of all X^2 values (i.e. sum(x1^2, x2^2, x3^2, ...))

Then when you need to compute the variance after 1,000th x values all you do is return

Var(X) = SumOfSquares/n - (Sum/n)^2

(this formula is given because Var(X) = E[X^2] - E[X]^2)

To compute Cov online

Keep track of these quantities (As mentioned above this means subtracting off quantities that are no longer in your window)

  1. Sum of all X values in your window (i.e. sum(x1, x2, x3, ....))
  2. Sum of all Y values in your window (i.e. sum(x1, x2, x3, ....))
  3. Sum of all X*Y values in your window (i..e sum(x1*y1, x2*y2, x3*y3))

All that is left is a little algebra...which I'll leave to you.

From a coder's point of view

I suggest writing a function like this:

double covar(cumulativeSumOfXs, cumulativeSumOfYs, cumulativeSumofXYs){
    //only simple algebra is necessary if these are your inputs
  • Thanks, so for sample covariance the formula I have is 1 / (n - 1) * [Sxy - (1 / n) * Sx * Sy]. My concern though is that Wikipedia warns this is subject to catastrophic cancellation. – MgSam Dec 6 '16 at 22:25
  • @MgSam -- Regarding catastrophic cancellation: many mathematical computations can "go wrong". The mere fact that a computation can go wrong doesn't necessarily mean you shouldn't use the "vulnerable" algorithm by default. For example, think of gaussian elimination (a simple method for solving linear systems of equations -- Linear Algebra 101 stuff). This extremely simple approach that works well much of the time. However, this simple approach works poorly when a matrix condition number is large. Should we throw-out a simple, easily understood algorithm just because bad inputs can exist? – Ivan Dec 7 '16 at 13:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.