How to efficiently calculate a moving Standard Deviation

Question

Below you can see my C# method to calculate Bollinger Bands for each point (moving average, up band, down band).

As you can see this method uses 2 for loops to calculate the moving standard deviation using the moving average. It used to contain an additional loop to calculate the moving average over the last n periods. This one I could remove by adding the new point value to total_average at the beginning of the loop and removing the i - n point value at the end of the loop.

My question now is basically: Can I remove the remaining inner loop in a similar way I managed with the moving average?

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

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

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

                for (int x = i; x > (i - period); x--)
                {
                    total_bollinger += Math.Pow(data.Values[x]["close"] - average, 2);
                }

                double stdev = Math.Sqrt(total_bollinger / 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"];
            }
        }
    }

Impatience

Answer 1

The answer is yes, you can. In the mid-80's I developed just such an algorithm (probably not original) in FORTRAN for a process monitoring and control application. Unfortunately, that was over 25 years ago and I do not remember the exact formulas, but the technique was an extension of the one for moving averages, with second order calculations instead of just linear ones.

---

After looking at your code some, I am think that I can suss out how I did it back then. Notice how your inner loop is making a Sum of Squares?:

            for (int x = i; x > (i - period); x--)
            {
                total_bollinger += Math.Pow(data.Values[x]["close"] - average, 2);
            }

in much the same way that your average must have originally had a Sum of Values? The only two differences are the order (its power 2 instead of 1) and that you are subtracting the average each value before you square it. Now that might look inseparable, but in fact they can be separated:

SUM(i=1; n){ (v[i] - k)^2 }

is

SUM(i=1..n){v[i]^2 -2*v[i]*k + k^2}

which becomes

SUM(i=1..n){v[i]^2 -2*v[i]*k} + k^2*n

which is

SUM(i=1..n){v[i]^2} + SUM(i=1..n){-2*v[i]*k} + k^2*n

which is also

SUM(i=1..n){v[i]^2} + SUM(i=1..n){-2*v[i]}*k + k^2*n

Now the first term is just a Sum of Squares, you handle that in the same way that you do the sum of Values for the average. The last term (k^2*n) is just the average squared times the period. Since you divide the result by the period anyway, you can just add the new average squared without the extra loop.

Finally, in the second term (SUM(-2*v[i]) * k), since SUM(v[i]) = total = k*n you can then change it into this:

-2 * k * k * n

or just -2*k^2*n, which is -2 times the average squared, once the period (n) is divided out again. So the final combined formula is:

SUM(i=1..n){v[i]^2} - n*k^2

or

SUM(i=1..n){values[i]^2} - period*(average^2)

(be sure to check the validity of this, since I am deriving it off the top of my head)

And incorporating into your code should look something like 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);
        }
    }
}

Decrescent

Thanks a lot! I was staring blind on this one. You only forgot to reduce the total_squares at the end: total_squares -= Math.Pow(data.Values[i - period + 1]["close"], 2); – Impatience

@odyth Nice reference! I had not realized this was in Knuth. I had actually read TAoCP several years before I wrote this in the 80's, and now I wonder if I subconsciously plagiarized it. – Decrescent

thanks for this, have been struggling to find a simple bollie band implementation that didn't have 2k lines of code :). this is brief and to the point. appreciated – Italy

Seems total_bollinger is never used? – Aboulia

Answer 2

The problem with approaches that calculate the sum of squares is that it and the square of sums can get quite large, and the calculation of their difference may introduce a very large error, so let's think of something better. For why this is needed, see the Wikipedia article on Algorithms for computing variance and John Cook on Theoretical explanation for numerical results)

First, instead of calculating the stddev let's focus on the variance. Once we have the variance, stddev is just the square root of the variance.

Suppose the data are in an array called x; rolling an n-sized window by one can be thought of as removing the value of x[0] and adding the value of x[n]. Let's denote the averages of x[0]..x[n-1] and x[1]..x[n] by µ and µ’ respectively. The difference between the variances of x[0]..x[n-1] and x[1]..x[n] is, after canceling out some terms and applying (a²-b²) = (a+b)(a-b):

Var[x[1],..,x[n]] - Var[x[0],..,x[n-1]] 
= (\sum_1^n x[i]² - n µ’²)/(n-1) - (\sum_0^{n-1} x[i]² - n µ²)/(n-1)
= (x[n]² - x[0]² - n(µ’² - µ²))/(n-1) 
= (x[n]-µ’ + x[0]-µ)(x[n]-x[0])/(n-1)

Therefore the variance is perturbed by something that doesn't require you to maintain the sum of squares, which is better for numerical accuracy.

You can calculate the mean and variance once in the beginning with a proper algorithm (Welford's method). After that, every time you have to replace a value in the window x[0] by another x[n] you update the average and variance like this:

new_Avg = Avg + (x[n]-x[0])/n
new_Var = Var + (x[n]-new_Avg + x[0]-Avg)(x[n] - x[0])/(n-1)
new_StdDev = sqrt(new_Var)

Manella

Thanks for this. I used it as the basis of an implementation in C# for the CLR. I discovered that, in practice, you may update such that new_Var is a very small negative number, and the sqrt fails. I introduced an if to limit the value to zero for this case. Not idea, but stable. This occurred when every value in my window had the same value (I used a window size of 20 and the value in question was 0.5, in case someone wants to try and reproduce this.) – Brion

Nice work! How would skewness or kurtosis would look like? Can it be derived the same way? – Aneurysm

FYI, i've implemented your algorithm in python (for the moment) and published it in this repo: github.com/ankostis/rollstats – Counts

Thank you for this answer! I used this answer to improve the performance of the Golang package for this answer: https://mcmap.net/q/54364/-peak-signal-detection-in-realtime-timeseries-data Golang package: github.com/MicahParks/peakdetect – Otherwise

Answer 3

Most important information has already been given above --- but maybe this is still of general interest.

A tiny Java library to calculate moving average and standard deviation is available here: https://github.com/tools4j/meanvar

The implementation is based on a variant of Welford's method mentioned above. Methods to remove and replace values have been derived that can be used for moving value windows.

Disclaimer: I am the author of the said library.

Payday

It's good form to disclose when associated with cited projects (if you are the co-contributor "Marco" of the project that is) – Putup

@Putup I couldn't agree more. I have updated the answer accordingly. – Payday

Answer 4

I've used commons-math (and contributed to that library!) for something very similar to this. It's open-source, porting to C# should be easy as store-bought pie (have you tried making a pie from scratch!?). Check it out: http://commons.apache.org/math/api-3.1.1/index.html. They have a StandardDeviation class. Go to town!

Southeastwardly

Thanks, but I don't think I need a math library for a simple calculation like this one. – Impatience

You're welcome! Sorry I didn't have the answer you're looking for. I definitely didn't mean to suggest porting the entire library! Just the minimum necessary code, which should be a few hundred lines or so. Note that I have no idea what legal / copyright restrictions apache has on that code, so you'd have to check that out. In case you pursue it, here is the link. So that + Variance + FastMath? – Southeastwardly

Answer 5

I just did it with Data From Binance Future API Hope this helps:

    using Newtonsoft.Json;
using System;
using System.Collections.Generic;
using System.Collections.ObjectModel;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.IO;
using System.Linq;
using System.Net;
using System.Text;
using System.Threading.Tasks;
using System.Windows.Forms;
using static System.Windows.Forms.VisualStyles.VisualStyleElement;

namespace Trading_Bot_1
{
    public class BOLL
    {
        private BollingerBandData graphdata = new BollingerBandData();
        private List<TickerData> data = new List<TickerData>();

        public BOLL(string url)
        {


            string js = getJsonFromUrl(url);
            //dynamic data = JObject.Parse(js);
            object[][] arrays = JsonConvert.DeserializeObject<object[][]>(js);

            data = new List<TickerData>();
            for (int i = 1; i < 400; i++)
            {
                data.Add(new TickerData
                {
                    Date = DateTime.Now,
                    Open = Convert.ToDouble(arrays[arrays.Length - i][1]),
                    High = Convert.ToDouble(arrays[arrays.Length - i][2]),
                    Low = Convert.ToDouble(arrays[arrays.Length - i][3]),
                    Close = Convert.ToDouble(arrays[arrays.Length - i][4]),
                    Volume = Math.Round(Convert.ToDouble(arrays[arrays.Length - i][4]), 0),
                    AdjClose = Convert.ToDouble(arrays[arrays.Length - i][6])
                });
            }

            graphdata.LowerBand.Add(1);
            graphdata.LowerBand.Add(2);
            graphdata.LowerBand.Add(3);
            graphdata.LowerBand.Add(1);

            graphdata.UpperBand.Add(1);
            graphdata.UpperBand.Add(2);
            graphdata.UpperBand.Add(3);
            graphdata.UpperBand.Add(4);

            graphdata.MovingAverageWindow.Add(10);
            graphdata.MovingAverageWindow.Add(20);
            graphdata.MovingAverageWindow.Add(40);
            graphdata.MovingAverageWindow.Add(50);

            graphdata.Length.Add(10);
            graphdata.Length.Add(30);
            graphdata.Length.Add(50);
            graphdata.Length.Add(100);

           // DataContext = graphdata;
        }

        public static string getJsonFromUrl(string url1)
        {
            var uri = String.Format(url1);
            WebClient client = new WebClient();
            client.UseDefaultCredentials = true;
            var data = client.DownloadString(uri);
            return data;
        }

        List<double> UpperBands = new List<double>();
        List<double> LowerBands = new List<double>();

        public List<List<double>> GetBOLLDATA(int decPlaces)
        {

            int datalength = graphdata.SelectedMovingAverage + graphdata.SelectedLength;

            string bands = "";
            for (int i = graphdata.SelectedLength - 1; i >= 0; i--)
            {
                List<double> price = new List<double>();
                for (int j = 0; j < graphdata.SelectedMovingAverage; j++)
                {
                    price.Add(data[i + j].Close);
                }
                double sma = CalculateAverage(price.ToArray());
                double sigma = CalculateSTDV(price.ToArray());
                double lower = sma - (graphdata.SelectedLowerBand * sigma);
                double upper = sma + (graphdata.SelectedUpperBand * sigma);
                UpperBands.Add(Math.Round( upper,decPlaces));
                LowerBands.Add(Math.Round(lower, decPlaces));
                bands += (Math.Round(upper, decPlaces) + " / " + Math.Round(lower, decPlaces)) + Environment.NewLine;
             //   graphdata.ChartData.Add(new ChartData() { SMA = sma, LowerBandData = lower, UpperBandData = upper });
            }
            //MessageBox.Show(bands);
            return new List<List<double>> { UpperBands, LowerBands };
        }
        public double[] GetBOLLDATA(int decPlaces, string a)
        {
            List<double> price = new List<double>();
            for (int j = 0; j < graphdata.SelectedMovingAverage; j++)
            {
                price.Add(data[j].Close);
            }
            double sma = CalculateAverage(price.ToArray());
            double sigma = CalculateSTDV(price.ToArray());
            double lower = sma - (graphdata.SelectedLowerBand * sigma);
            double upper = sma + (graphdata.SelectedUpperBand * sigma);
            return new double[] { Math.Round(upper, decPlaces), Math.Round(lower, decPlaces) };
        }


        private double CalculateAverage(double[] data)
        {
            int count = data.Length;
            double sum = 0;

            for (int i = 0; i < count; i++)
            {
                sum += data[i];
            }

            return sum / count;
        }

        private double CalculateVariance(double[] data)
        {
            int count = data.Length;
            double sum = 0;
            double avg = CalculateAverage(data);

            for (int i = 0; i < count; i++)
            {
                sum += (data[i] - avg) * (data[i] - avg);
            }

            return sum / (count - 1);
        }

        private double CalculateSTDV(double[] data)
        {
            double var = CalculateVariance(data);

            return Math.Sqrt(var);
        }

    }
    public class ChartData
    {
        public double UpperBandData
        { get; set; }

        public double LowerBandData
        { get; set; }

        public double SMA
        { get; set; }
    }

    public class BollingerBandData : INotifyPropertyChanged
    {
        private ObservableCollection<int> _lowerBand;
        private ObservableCollection<int> _upperBand;
        private ObservableCollection<int> _movingAvg;
        private ObservableCollection<int> _length;
        private ObservableCollection<ChartData> _chartData;
        private int _selectedLowerBand;
        private int _selectedUpperBand;
        private int _selectedMovingAvg;
        private int _selectedLength;

        public BollingerBandData()
        {
            _lowerBand = new ObservableCollection<int>();
            _upperBand = new ObservableCollection<int>();
            _movingAvg = new ObservableCollection<int>();
            _length = new ObservableCollection<int>();
            _chartData = new ObservableCollection<ChartData>();

            SelectedLowerBand = 2;
            SelectedUpperBand = 2;
            SelectedMovingAverage = 20;
            SelectedLength = 5;
        }

        public ObservableCollection<ChartData> ChartData
        {
            get
            {
                return _chartData;
            }
            set
            {
                _chartData = value;
                RaisePropertyChanged("ChartData");
            }
        }

        public ObservableCollection<int> LowerBand
        {
            get
            {
                return _lowerBand;
            }
            set
            {
                _lowerBand = value;
                RaisePropertyChanged("LowerBand");
            }
        }

        public ObservableCollection<int> UpperBand
        {
            get
            {
                return _upperBand;
            }
            set
            {
                _upperBand = value;
                RaisePropertyChanged("UpperBand");
            }
        }

        public ObservableCollection<int> MovingAverageWindow
        {
            get
            {
                return _movingAvg;
            }
            set
            {
                _movingAvg = value;
                RaisePropertyChanged("MovingAverageWindow");
            }
        }

        public ObservableCollection<int> Length
        {
            get
            {
                return _length;
            }
            set
            {
                _length = value;
                RaisePropertyChanged("Length");
            }
        }
        public int SelectedLowerBand
        {
            get
            {
                return _selectedLowerBand;
            }
            set
            {
                _selectedLowerBand = value;
                RaisePropertyChanged("SelectedLowerBand");
            }
        }

        public int SelectedUpperBand
        {
            get
            {
                return _selectedUpperBand;
            }
            set
            {
                _selectedUpperBand = value;
                RaisePropertyChanged("SelectedUpperBand");
            }
        }

        public int SelectedMovingAverage
        {
            get
            {
                return _selectedMovingAvg;
            }
            set
            {
                _selectedMovingAvg = value;
                RaisePropertyChanged("SelectedMovingAverage");
            }
        }

        public int SelectedLength
        {
            get
            {
                return _selectedLength;
            }
            set
            {
                _selectedLength = value;
                RaisePropertyChanged("SelectedLength");
            }
        }

        public event PropertyChangedEventHandler PropertyChanged;

        private void RaisePropertyChanged(string propertyName)
        {
            PropertyChangedEventHandler handler = this.PropertyChanged;

            if (handler != null)
            {
                handler(this, new PropertyChangedEventArgs(propertyName));
            }
        }
    }

    public class TickerData
    {
        public DateTime Date
        { get; set; }

        public double Open
        { get; set; }

        public double High
        { get; set; }

        public double Low
        { get; set; }

        public double Close
        { get; set; }

        public double Volume
        { get; set; }

        public double AdjClose
        { get; set; }
    }




}

Caterwaul