Calculate a moving Standard Deviation

Question

i found the following code snippet here on stackoverflow, but i have the problem that the stdev becomes NaN. Any ideas how to fix this?

public static void AddBollingerBands(ref SortedList<DateTime, Dictionary<string, double>> data, int period, int factor)
{
    double total_average = 0;
    double total_squares = 0;

    for (int i = 0; i < data.Count(); i++)
    {
        total_average += data.Values[i]["close"];
        total_squares += Math.Pow(data.Values[i]["close"], 2);

        if (i >= period - 1)
        {
            double total_bollinger = 0;
            double average = total_average / period;

            double stdev = Math.Sqrt((total_squares - Math.Pow(total_average,2)/period) / period);
            data.Values[i]["bollinger_average"] = average;
            data.Values[i]["bollinger_top"] = average + factor * stdev;
            data.Values[i]["bollinger_bottom"] = average - factor * stdev;

            total_average -= data.Values[i - period + 1]["close"];
            total_squares -= Math.Pow(data.Values[i - period + 1]["close"], 2);
        }
    }
}

Answer 1

A numerically more stable variant is preferable when doing incremental / moving average and standard deviation calculations. One way to do this is using Knuth's algorithm, as shown in the code block below:

public class MovingAverageCalculator
{
    public MovingAverageCalculator(int period)
    {
        _period = period;
        _window = new double[period];
    }

    public double Average 
    {
        get { return _average; }
    }

    public double StandardDeviation
    {
        get 
        {
            var variance = Variance;
            if (variance >= double.Epsilon)
            {
                var sd = Math.Sqrt(variance);
                return double.IsNaN(sd) ? 0.0 : sd;
            }
            return 0.0;
        }
    }

    public double Variance
    {
        get 
        { 
            var n = N;
            return n > 1 ? _variance_sum / (n - 1) : 0.0; 
        }
    }

    public bool HasFullPeriod
    {
        get { return _num_added >= _period; }
    }

    public IEnumerable<double> Observations
    {
        get { return _window.Take(N); }
    }

    public int N
    {
        get { return Math.Min(_num_added, _period); }
    }

    public void AddObservation(double observation)
    {
        // Window is treated as a circular buffer.
        var ndx = _num_added % _period;
        var old = _window[ndx];     // get value to remove from window
        _window[ndx] = observation; // add new observation in its place.
        _num_added++;

        // Update average and standard deviation using deltas
        var old_avg = _average;
        if (_num_added <= _period)
        {
            var delta = observation - old_avg;
            _average += delta / _num_added;
            _variance_sum += (delta * (observation - _average));
        } 
        else // use delta vs removed observation.
        {
            var delta = observation - old;
            _average += delta / _period;
            _variance_sum += (delta * ((observation - _average) + (old - old_avg)));
        }
    }

    private readonly int _period;
    private readonly double[] _window;
    private int _num_added;
    private double _average;
    private double _variance_sum;
}

You could then use it in the following manner in your code example:

public static void AddBollingerBands(ref SortedList<DateTime, Dictionary<string, double>> data, int period, int factor)
{
    var moving_avg = new MovingAverageCalculator(period);
    for (int i = 0; i < data.Count(); i++)
    {
        moving_avg.AddObservation(data.Values[i]["close"]);
        if (moving_avg.HasFullPeriod)
        {
            var average = moving_avg.Average;
            var limit = factor * moving_avg.StandardDeviation;
            data.Values[i]["bollinger_average"] = average;
            data.Values[i]["bollinger_top"] = average + limit;
            data.Values[i]["bollinger_bottom"] = average - limit;
        }
    }
}