rolling.hpp

Go to the documentation of this file.

#pragma once
 
#include "./numeric.hpp"
 
#include <cmath>
#include <cstddef>
#include <limits>
#include <stdexcept>
#include <string>
 
namespace algotradelib {
namespace analytics {

 
namespace rolling_detail {

struct Reduction {
  std::size_t count_{0};
  Scalar sum_{0};
};

inline Reduction reduce(const math::Observed& parts, math::Index first, math::Index length) {
  Reduction result;
  result.count_ = static_cast<std::size_t>(parts.present_.segment(first, length).sum());
  result.sum_ = parts.values_.segment(first, length).sum();
  return result;
}

struct Dispersion {
  bool reported_{false};
  std::size_t count_{0};
  Scalar mean_{0};
  Scalar deviation_{0};
};

inline Dispersion dispersion(
  const math::Observed& parts,
  math::Index first,
  math::Index length,
  std::size_t needed
) {
  Dispersion result;
  const auto reduced = reduce(parts, first, length);
  result.count_ = reduced.count_;
  if (reduced.count_ < needed || reduced.count_ < 2) {
    return result;
  }
  result.mean_ = reduced.sum_ / static_cast<Scalar>(reduced.count_);
  const Scalar squared = (parts.values_.segment(first, length) -
    result.mean_ * parts.present_.segment(first, length)).squaredNorm();
  result.deviation_ = std::sqrt(squared / static_cast<Scalar>(reduced.count_ - 1));
  result.reported_ = true;
  return result;
}

inline void check(std::size_t window, std::size_t min_observed) {
  if (window == 0) {
    throw std::invalid_argument("rolling window must be at least one row");
  }
  if (min_observed == 0) {
    throw std::invalid_argument("a window must require at least one observation");
  }
  if (min_observed > window) {
    throw std::invalid_argument(
      "rolling window of " + std::to_string(window) + " rows cannot hold " +
      std::to_string(min_observed) + " observations"
    );
  }
}

inline math::Index window_start(math::Index row, std::size_t window) {
  const auto span = static_cast<math::Index>(window);
  return row + 1 >= span ? row + 1 - span : math::Index{0};
}
 
} // namespace rolling_detail

inline math::Vector rolling_sum(
  math::SeriesView values,
  std::size_t window,
  std::size_t min_observed
) {
  rolling_detail::check(window, min_observed);
  const auto parts = math::observed_parts(values);
  math::Vector result = math::unobserved(static_cast<std::size_t>(values.size()));
  for (math::Index row = 0; row < values.size(); ++row) {
    const math::Index first = rolling_detail::window_start(row, window);
    const auto reduced = rolling_detail::reduce(parts, first, row - first + 1);
    if (reduced.count_ >= min_observed) {
      result[row] = reduced.sum_;
    }
  }
  return result;
}

inline math::Vector rolling_sum(math::SeriesView values, std::size_t window) {
  return rolling_sum(values, window, window);
}

inline math::Vector rolling_mean(
  math::SeriesView values,
  std::size_t window,
  std::size_t min_observed
) {
  rolling_detail::check(window, min_observed);
  const auto parts = math::observed_parts(values);
  math::Vector result = math::unobserved(static_cast<std::size_t>(values.size()));
  for (math::Index row = 0; row < values.size(); ++row) {
    const math::Index first = rolling_detail::window_start(row, window);
    const auto reduced = rolling_detail::reduce(parts, first, row - first + 1);
    if (reduced.count_ >= min_observed) {
      result[row] = reduced.sum_ / static_cast<Scalar>(reduced.count_);
    }
  }
  return result;
}

inline math::Vector rolling_mean(math::SeriesView values, std::size_t window) {
  return rolling_mean(values, window, window);
}

inline math::Vector rolling_stdev(
  math::SeriesView values,
  std::size_t window,
  std::size_t min_observed
) {
  rolling_detail::check(window, min_observed);
  const auto parts = math::observed_parts(values);
  math::Vector result = math::unobserved(static_cast<std::size_t>(values.size()));
  for (math::Index row = 0; row < values.size(); ++row) {
    const math::Index first = rolling_detail::window_start(row, window);
    const auto spread =
      rolling_detail::dispersion(parts, first, row - first + 1, min_observed);
    if (spread.reported_) {
      result[row] = spread.deviation_;
    }
  }
  return result;
}

inline math::Vector rolling_stdev(math::SeriesView values, std::size_t window) {
  return rolling_stdev(values, window, window);
}

inline math::Vector rolling_min(
  math::SeriesView values,
  std::size_t window,
  std::size_t min_observed
) {
  rolling_detail::check(window, min_observed);
  const auto parts = math::observed_parts(values);
  // Absent rows become `+inf` so they lose every comparison without needing a mask inside the
  // window. A window with no observation at all is rejected by the count, never read from here.
  const math::Vector ceiling =
    math::observed_or(values, std::numeric_limits<Scalar>::infinity());
  math::Vector result = math::unobserved(static_cast<std::size_t>(values.size()));
  for (math::Index row = 0; row < values.size(); ++row) {
    const math::Index first = rolling_detail::window_start(row, window);
    const math::Index length = row - first + 1;
    if (parts.present_.segment(first, length).sum() >= static_cast<Scalar>(min_observed)) {
      result[row] = ceiling.segment(first, length).minCoeff();
    }
  }
  return result;
}

inline math::Vector rolling_min(math::SeriesView values, std::size_t window) {
  return rolling_min(values, window, window);
}

inline math::Vector rolling_max(
  math::SeriesView values,
  std::size_t window,
  std::size_t min_observed
) {
  rolling_detail::check(window, min_observed);
  const auto parts = math::observed_parts(values);
  const math::Vector floor =
    math::observed_or(values, -std::numeric_limits<Scalar>::infinity());
  math::Vector result = math::unobserved(static_cast<std::size_t>(values.size()));
  for (math::Index row = 0; row < values.size(); ++row) {
    const math::Index first = rolling_detail::window_start(row, window);
    const math::Index length = row - first + 1;
    if (parts.present_.segment(first, length).sum() >= static_cast<Scalar>(min_observed)) {
      result[row] = floor.segment(first, length).maxCoeff();
    }
  }
  return result;
}

inline math::Vector rolling_max(math::SeriesView values, std::size_t window) {
  return rolling_max(values, window, window);
}

inline math::Vector sma(math::SeriesView values, std::size_t window, std::size_t min_observed) {
  return rolling_mean(values, window, min_observed);
}

inline math::Vector sma(math::SeriesView values, std::size_t window) {
  return rolling_mean(values, window, window);
}

inline math::Vector wma(
  math::SeriesView values,
  std::size_t window,
  std::size_t min_observed
) {
  rolling_detail::check(window, min_observed);
  const auto parts = math::observed_parts(values);
  math::Vector weights(static_cast<math::Index>(window));
  for (math::Index at = 0; at < weights.size(); ++at) {
    weights[at] = static_cast<Scalar>(at + 1);
  }
 
  math::Vector result = math::unobserved(static_cast<std::size_t>(values.size()));
  for (math::Index row = 0; row < values.size(); ++row) {
    const math::Index first = rolling_detail::window_start(row, window);
    const math::Index length = row - first + 1;
    if (parts.present_.segment(first, length).sum() < static_cast<Scalar>(min_observed)) {
      continue;
    }
    const auto used = weights.head(length);
    const Scalar divisor = parts.present_.segment(first, length).dot(used);
    if (divisor > Scalar{0}) {
      result[row] = parts.values_.segment(first, length).dot(used) / divisor;
    }
  }
  return result;
}

inline math::Vector wma(math::SeriesView values, std::size_t window) {
  return wma(values, window, window);
}

inline Scalar alpha_from_span(std::size_t span) {
  if (span == 0) {
    throw std::invalid_argument("an ema span must be at least one row");
  }
  return Scalar{2} / (static_cast<Scalar>(span) + Scalar{1});
}
 
namespace rolling_detail {

inline void check_alpha(Scalar alpha) {
  if (!(Scalar{0} < alpha && alpha <= Scalar{1})) {
    throw std::invalid_argument("an ema smoothing factor must be in (0, 1]");
  }
}
 
} // namespace rolling_detail

inline Scalar ema_step(Scalar level, Scalar value, Scalar alpha) {
  rolling_detail::check_alpha(alpha);
  if (!observed(value)) {
    return level;
  }
  return observed(level) ? level + alpha * (value - level) : value;
}

Scalar ema_step(Scalar level, Scalar value, std::size_t span) = delete;

inline math::Vector ema(math::SeriesView values, Scalar alpha, std::size_t min_observed = 1) {
  rolling_detail::check_alpha(alpha);
  if (min_observed == 0) {
    throw std::invalid_argument("ema must require at least one observation");
  }
  math::Vector result = math::unobserved(static_cast<std::size_t>(values.size()));
  Scalar level{0};
  std::size_t count = 0;
  for (math::Index row = 0; row < values.size(); ++row) {
    const Scalar value = values[row];
    if (!observed(value)) {
      continue;
    }
    // Seeded by count rather than by an unobserved level, so a level an infinity turned into `NaN`
    // stays `NaN` instead of being mistaken for "no level yet".
    level = count == 0 ? value : level + alpha * (value - level);
    ++count;
    if (count >= min_observed) {
      result[row] = level;
    }
  }
  return result;
}

math::Vector ema(math::SeriesView values, std::size_t span, std::size_t min_observed = 1) = delete;
 
} // namespace analytics
} // namespace algotradelib