Incremental, weighted average, variance, covariance

incremental.test.ts

import { round } from "mathjs"
import { average, sampleCovariance } from "simple-statistics"

const sum = (array: number[]) => {
  return array.reduce((a, b) => a + b, 0)
}

const weightAverage = (prices: number[], weights: number[]) => {
  const wSum = sum(weights)
  const weightedPrices = prices.map((price, index) => {
    return price * weights[index]
  })
  return sum(weightedPrices) / wSum
}

const incrementalWeightedAverage = (
  prevAverage: number,
  prevWeightSum: number,
  newPrice: number,
  newWeight: number,
) => {
  const newWeightSum = prevWeightSum + newWeight
  const newAverage =
    prevAverage + ((newPrice - prevAverage) * newWeight) / newWeightSum
  return { newAverage, newWeightSum }
}

describe("average test", () => {
  test("incrementalWeightedAverage", () => {
    const x = [6, 2, 3, 4, 5]
    const weights = [0.5, 1, 0.5, 1, 0.5]
    let prevAverage = 0
    let prevWeightSum = 0
    for (let i = 0; i < x.length; i++) {
      const result = incrementalWeightedAverage(
        prevAverage,
        prevWeightSum,
        x[i],
        weights[i],
      )
      prevAverage = result.newAverage
      prevWeightSum = result.newWeightSum

      const sliced = x.slice(0, i + 1)
      const slicedWeights = weights.slice(0, i + 1)
      const correctNewWeightedAverage = weightAverage(sliced, slicedWeights)
      // eslint-disable-next-line
      if (false)
        console.log({
          sliced,
          slicedWeights,
          correctNewWeightedAverage,
          result,
        })
      expect(result.newAverage).toBe(correctNewWeightedAverage)
    }
  })
})

const weightStandardDeviation = (prices: number[], weights: number[]) => {
  const covariance = weightVariance(prices, weights)
  return Math.sqrt(covariance)
}

export const weightVariance = (prices: number[], weights: number[]) => {
  const wSum = sum(weights)
  const wAverage = weightAverage(prices, weights)
  const weightedPrices = prices.map((price, index) => {
    return weights[index] * (price - wAverage) ** 2
  })
  return sum(weightedPrices) / wSum
}

/**
 * https://fanf2.user.srcf.net/hermes/doc/antiforgery/stats.pdf
 */
const incrementalWeightedVariance = (
  prevVariance: number,
  prevAverage: number,
  newAverage: number,
  newPrice: number,
  newWeight: number,
) => {
  const newVariance =
    prevVariance +
    newWeight * (newPrice - prevAverage) * (newPrice - newAverage)
  return { newVariance }
}

describe("variance test", () => {
  test("weightVariance", () => {
    const x = [6, 2, 3, 4, 5]
    const weightsA = [0.5, 1, 0.5, 1, 0.5]
    const weightsB = [1, 2, 1, 2, 1]
    expect(weightVariance(x, weightsA)).toBe(weightVariance(x, weightsB))
  })

  test("incrementalWeightedVariance", () => {
    const x = [6, 2, 3, 4, 5]
    const weights = [0.5, 1, 0.5, 1, 0.5]
    let prevAverage = 0
    let prevVariance = 0
    let prevWeightSum = 0
    for (let i = 0; i < x.length; i++) {
      const averageResult = incrementalWeightedAverage(
        prevAverage,
        prevWeightSum,
        x[i],
        weights[i],
      )
      const newPrice = x[i]
      const newWeight = weights[i]
      const result = incrementalWeightedVariance(
        prevVariance,
        prevAverage,
        averageResult.newAverage,
        newPrice,
        newWeight,
      )
      prevAverage = averageResult.newAverage
      prevWeightSum = averageResult.newWeightSum
      prevVariance = result.newVariance

      const standardDeviation = Math.sqrt(
        result.newVariance / averageResult.newWeightSum,
      )

      const sliced = x.slice(0, i + 1)
      const slicedWeights = weights.slice(0, i + 1)
      // eslint-disable-next-line
      if (false)
        console.log({
          standardDeviation: round(standardDeviation, 15),
          sliced,
          prevAverage,
          prevVariance,
          newPrice,
          newWeight,
          result,
        })
      if (i === 0) continue
      expect(round(standardDeviation, 15)).toBe(
        round(weightStandardDeviation(sliced, slicedWeights), 15),
      )
    }
  })
})

/**
 * https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Weighted_batched_version
 */
const incrementalWeightedCovariance = ({
  prevCovarianceSum,
  prevWeightSum,
  prevAverageY,
  newAverageX,
  newX,
  newY,
  newWeight,
}: {
  prevCovarianceSum: number
  prevWeightSum: number
  prevAverageY: number
  newAverageX: number
  newX: number
  newY: number
  newWeight: number
}) => {
  const newWeightSum = prevWeightSum + newWeight
  const newCovarianceSum =
    prevCovarianceSum + newWeight * (newX - newAverageX) * (newY - prevAverageY)
  const newCovariance = newCovarianceSum / prevWeightSum // why prevWeightSum works instead of newWeightSum?
  return { newCovariance, newWeightSum, newCovarianceSum }
}

describe("covariance test", () => {
  test("incrementalWeightedCovariance", () => {
    const data = Array.from({ length: 20 }, () => [
      Math.random() * 100,
      Math.random() * 100,
    ])
    const x = data.map((d) => d[0])
    const y = data.map((d) => d[1])
    const weights = Array.from({ length: data.length }, () => 1) // TODO: test with weights

    let prevCovarianceSum = 0
    let prevAverageY = 0
    let prevWeightSum = 0

    for (let i = 0; i < data.length; i++) {
      const slicedX = x.slice(0, i + 1)
      const slicedY = y.slice(0, i + 1)
      const slicedWeights = weights.slice(0, i + 1)
      const newAverageX = average(slicedX)
      const result = incrementalWeightedCovariance({
        prevCovarianceSum,
        prevWeightSum,
        prevAverageY,
        newAverageX,
        newX: x[i],
        newY: y[i],
        newWeight: weights[i],
      })
      // eslint-disable-next-line
      if (false)
        console.log({
          prevCovarianceSum,
          prevWeightSum,
          slicedX,
          slicedY,
          slicedWeights,
          prevAverageY,
          newAverageX,
          newValues: [x[i], y[i], weights[i]],
          result,
        })
      if (i > 0) {
        const correctCovariance = sampleCovariance(slicedX, slicedY)
        expect(round(result.newCovariance, 10)).toBe(
          round(correctCovariance, 10),
        )
      }

      prevCovarianceSum = result.newCovarianceSum
      prevAverageY = average(slicedY)
      prevWeightSum = result.newWeightSum
    }
  })
})