marketdata/MarketDataLib/Numerics/LeastSquares.cs

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace MarketData.Numerical
{
  public class LeastSquaresResult
  {
    public double[] LeastSquares { get; set; }
    public double YIntercept { get; set; }
    public double XIntercept { get; set; }
    public double Slope { get; set; }
    public void Extend(int itemCount, double slope)      // extend the result set by itemCount by applying the slope to subsequent values
    {
      if (null == LeastSquares || 0 == LeastSquares.Length) return;
      double[] extendedItems = new double[itemCount + LeastSquares.Length];
      Array.Copy(LeastSquares, 0, extendedItems, itemCount, LeastSquares.Length);
      for (int index = itemCount - 1; index >= 0; index--)
      {
        extendedItems[index] = extendedItems[index + 1] + slope;
      }
      LeastSquares = extendedItems;
    }
  }
  public class LeastSquaresResultWithR2 : LeastSquaresResult
  {
    public double RSquared { get; set; }
  }
// **********************************************************************************************************************************************
  public class LeastSquaresHelper
  {
    private LeastSquaresHelper()
    {
    }
// This version reverses the slope because the input data is in reverse order
// The input data should be in the following manner observations[count-1]=earliest data point, observations[0]=most recent data point
// This version should be retained as-is because it is being referenced in models etc.,
    public static LeastSquaresResult CalculateLeastSquares(double[] observations)
    {
      try
      {
        LeastSquaresResult leastSquaresResult = new LeastSquaresResult();
        double[] xSeries = new double[observations.Length];
        double[] xMinusXMean = new double[observations.Length];
        double[] yMinusYMean = new double[observations.Length];
        double[] meanProduct = new double[observations.Length];
        double[] xMinusXMeanSquared = new double[observations.Length];
        double[] leastSquares = new double[observations.Length];
        double xMean = double.NaN;
        double yMean = double.NaN;
        double slope = double.NaN;
        double yIntercept = double.NaN;
        for (int index = 0; index < xSeries.Length; index++) xSeries[index] = index + 1;
        xMean = Numerics.Mean(ref xSeries);
        yMean = Numerics.Mean(ref observations);
        for (int index = 0; index < observations.Length; index++)
        {
          xMinusXMean[index] = xSeries[index] - xMean;
          yMinusYMean[index] = observations[index] - yMean;
          meanProduct[index] = xMinusXMean[index] * yMinusYMean[index];
          xMinusXMeanSquared[index] = xMinusXMean[index] * xMinusXMean[index];
        }
        slope = Numerics.Sum(ref meanProduct) / Numerics.Sum(ref xMinusXMeanSquared);
        yIntercept = Numerics.YIntercept(xMean, yMean, slope);
        for (int index = 0; index < leastSquares.Length; index++) leastSquares[index] = slope * xSeries[index] + yIntercept;
        leastSquaresResult.LeastSquares = leastSquares;
        leastSquaresResult.YIntercept = yIntercept;
        leastSquaresResult.XIntercept = xMean;
        leastSquaresResult.Slope = slope * -1.00;
        return leastSquaresResult;
      }
      catch (Exception exception)
      {
        MDTrace.WriteLine(LogLevel.DEBUG, exception.ToString());
        return null;
      }
    }
// This version uses the true slope whereas the version above reverses the slope due to the arrangement of the input data
// The input data should be in the following manner observations[0]=earliest data point, observations[count-1]=most recent data point
    public static LeastSquaresResultWithR2 CalculateLeastSquaresWithR2(double[] observations)
    {
      try
      {
        LeastSquaresResultWithR2 leastSquaresResultWithR2 = new LeastSquaresResultWithR2();
        double[] xSeries = new double[observations.Length];
        double[] xMinusXMean = new double[observations.Length];
        double[] yMinusYMean = new double[observations.Length];
        double[] meanProduct = new double[observations.Length];
        double[] xMinusXMeanSquared = new double[observations.Length];
        double[] leastSquares = new double[observations.Length];
        double xMean = double.NaN;
        double yMean = double.NaN;
        double slope = double.NaN;
        double yIntercept = double.NaN;
        for (int index = 0; index < xSeries.Length; index++) xSeries[index] = index + 1;
        xMean = Numerics.Mean(ref xSeries);
        yMean = Numerics.Mean(ref observations);
        for (int index = 0; index < observations.Length; index++)
        {
          xMinusXMean[index] = xSeries[index] - xMean;
          yMinusYMean[index] = observations[index] - yMean;
          meanProduct[index] = xMinusXMean[index] * yMinusYMean[index];
          xMinusXMeanSquared[index] = xMinusXMean[index] * xMinusXMean[index];
        }
        slope = Numerics.Sum(ref meanProduct) / Numerics.Sum(ref xMinusXMeanSquared);
        yIntercept = Numerics.YIntercept(xMean, yMean, slope);
        for (int index = 0; index < leastSquares.Length; index++) leastSquares[index] = slope * xSeries[index] + yIntercept;
        leastSquaresResultWithR2.LeastSquares = leastSquares;
        leastSquaresResultWithR2.YIntercept = yIntercept;
        leastSquaresResultWithR2.XIntercept = xMean;
        leastSquaresResultWithR2.Slope = slope;
        // calculate the R2
        double[] estimated = new double[observations.Length];
        double[] estimatedLessMeanSquared = new double[observations.Length];
        double sumOfEstimatedLessMeanSquared = 0.00;
        double sumofsquares = 0.00;
        foreach (double observation in observations) sumofsquares += Math.Pow(observation - yMean, 2);
        for (int index = 0; index < observations.Length; index++) estimated[index] = yIntercept + (slope * (index + 1));
        for (int index = 0; index < observations.Length; index++) estimatedLessMeanSquared[index] = Math.Pow(estimated[index] - yMean, 2);
        foreach (double value in estimatedLessMeanSquared) sumOfEstimatedLessMeanSquared += value;
        leastSquaresResultWithR2.RSquared = sumOfEstimatedLessMeanSquared / sumofsquares;
        return leastSquaresResultWithR2;
      }
      catch (Exception exception)
      {
        MDTrace.WriteLine(LogLevel.DEBUG, exception.ToString());
        return null;
      }
    }
  }
}