Rolling variance algorithm

I'm trying to find an efficient, numerically stable algorithm to calculate a rolling variance (for instance, a variance over a 20-period rolling window). I'm aware of the Welford algorithm that efficiently computes the running variance for a stream of numbers (it requires only one pass), but am not sure if this can be adapted for a rolling window. I would also like the solution to avoid the accuracy problems discussed at the top of this article by John D. Cook. A solution in any language is fine.

Solution 1:

I've run across this problem as well. There are some great posts out there in computing the running cumulative variance such as John Cooke's Accurately computing running variance post and the post from Digital explorations, Python code for computing sample and population variances, covariance and correlation coefficient. Just could not find any that were adapted to a rolling window.

The Running Standard Deviations post by Subluminal Messages was critical in getting the rolling window formula to work. Jim takes the power sum of the squared differences of the values versus Welford’s approach of using the sum of the squared differences of the mean. Formula as follows:

PSA today = PSA(yesterday) + (((x today * x today) - x yesterday)) / n

• x = value in your time series
• n = number of values you've analyzed so far.

But, to convert the Power Sum Average formula to a windowed variety you need tweak the formula to the following:

PSA today = PSA yesterday + (((x today * x today) - (x yesterday * x Yesterday) / n

• x = value in your time series
• n = number of values you've analyzed so far.

You'll also need the Rolling Simple Moving Average formula:

SMA today = SMA yesterday + ((x today - x today - n) / n

• x = value in your time series
• n = period used for your rolling window.

From there you can compute the Rolling Population Variance:

Population Var today = (PSA today * n - n * SMA today * SMA today) / n

Or the Rolling Sample Variance:

Sample Var today = (PSA today * n - n * SMA today * SMA today) / (n - 1)

I've covered this topic along with sample Python code in a blog post a few years back, Running Variance.

Hope this helps.

Please note: I provided links to all the blog posts and math formulas in Latex (images) for this answer. But, due to my low reputation (< 10); I'm limited to only 2 hyperlinks and absolutely no images. Sorry about this. Hope this doesn't take away from the content.

Solution 2:

I have been dealing with the same issue.

Mean is simple to compute iteratively, but you need to keep the complete history of values in a circular buffer.

next_index = (index + 1) % window_size;    // oldest x value is at next_index, wrapping if necessary.

new_mean = mean + (x_new - xs[next_index])/window_size;

I have adapted Welford's algorithm and it works for all the values that I have tested with.

varSum = var_sum + (x_new - mean) * (x_new - new_mean) - (xs[next_index] - mean) * (xs[next_index] - new_mean);

xs[next_index] = x_new;
index = next_index;

To get the current variance just divide varSum by the window size: variance = varSum / window_size;

Solution 3:

If you prefer code over words (heavily based on DanS' post): http://calcandstuff.blogspot.se/2014/02/rolling-variance-calculation.html

public IEnumerable RollingSampleVariance(IEnumerable data, int sampleSize)
{
    double mean = 0;
    double accVar = 0;

    int n = 0;
    var queue = new Queue(sampleSize);

    foreach(var observation in data)
    {
        queue.Enqueue(observation);
        if (n < sampleSize)
        {
            // Calculating first variance
            n++;
            double delta = observation - mean;
            mean += delta / n;
            accVar += delta * (observation - mean);
        }
        else
        {
            // Adjusting variance
            double then = queue.Dequeue();
            double prevMean = mean;
            mean += (observation - then) / sampleSize;
            accVar += (observation - prevMean) * (observation - mean) - (then - prevMean) * (then - mean);
        }

        if (n == sampleSize)
            yield return accVar / (sampleSize - 1);
    }
}