indicators.hpp

Go to the documentation of this file.

#pragma once
 
#include "./rolling.hpp"
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <stdexcept>
 
namespace algotradelib {
namespace analytics {

 
namespace indicator_detail {

inline math::Vector wilder(math::SeriesView values, std::size_t period) {
  const auto span = static_cast<Scalar>(period);
  math::Vector result = math::unobserved(static_cast<std::size_t>(values.size()));
  Scalar average{0};
  Scalar seed{0};
  std::size_t count = 0;
  for (math::Index row = 0; row < values.size(); ++row) {
    const Scalar value = values[row];
    if (!observed(value)) {
      continue;
    }
    ++count;
    if (count < period) {
      seed += value;
      continue;
    }
    if (count == period) {
      average = (seed + value) / span;
    } else {
      average = (average * (span - Scalar{1}) + value) / span;
    }
    result[row] = average;
  }
  return result;
}

struct Paired {
  math::Vector left_;
  math::Vector right_;
  math::Vector present_;
};

inline Paired pair_series(math::SeriesView left, math::SeriesView right) {
  const auto left_parts = math::observed_parts(left);
  const auto right_parts = math::observed_parts(right);
  Paired paired;
  paired.present_ = math::both_present(left_parts.present_, right_parts.present_);
  paired.left_ = left_parts.values_.cwiseProduct(paired.present_);
  paired.right_ = right_parts.values_.cwiseProduct(paired.present_);
  return paired;
}

struct Pairing {
  bool reported_{false};
  std::size_t count_{0};
  Scalar mean_left_{0};
  Scalar mean_right_{0};
};

inline Pairing reduce_pair(
  const Paired& paired,
  math::Index first,
  math::Index length,
  std::size_t needed
) {
  Pairing result;
  const Scalar count = paired.present_.segment(first, length).sum();
  result.count_ = static_cast<std::size_t>(count);
  if (result.count_ < needed || result.count_ < 2) {
    return result;
  }
  result.mean_left_ = paired.left_.segment(first, length).sum() / count;
  result.mean_right_ = paired.right_.segment(first, length).sum() / count;
  result.reported_ = true;
  return result;
}

inline math::Vector centred(
  const math::Vector& values,
  const math::Vector& present,
  math::Index first,
  math::Index length,
  Scalar mean
) {
  return values.segment(first, length) - mean * present.segment(first, length);
}

inline void require_same_rows(math::SeriesView left, math::SeriesView right) {
  if (left.size() != right.size()) {
    throw std::invalid_argument("paired series must have the same number of rows");
  }
}
 
} // namespace indicator_detail

inline math::Vector zscore(
  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) {
    if (!observed(values[row])) {
      continue;
    }
    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_) {
      continue;
    }
    if (spread.deviation_ == Scalar{0} ||
        constant_detail::rounded_constant(
          spread.deviation_, spread.mean_, spread.count_, parts.values_, parts.present_, first, row - first + 1
        )) {
      // The data was there and the statistic does not exist for it, which is a different
      // answer from "no data yet" and now says so. Equal values that are not binary fractions
      // leave a deviation of a few ulps rather than zero; dividing by it would not be a z-score.
      result[row] = math::undefined();
      continue;
    }
    result[row] = (values[row] - spread.mean_) / spread.deviation_;
  }
  return result;
}

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

struct BollingerBands {
  math::Vector middle_;
  math::Vector upper_;
  math::Vector lower_;
};

inline BollingerBands bollinger(
  math::SeriesView values,
  std::size_t window,
  Scalar multiplier,
  std::size_t min_observed
) {
  rolling_detail::check(window, min_observed);
  const auto parts = math::observed_parts(values);
  const auto rows = static_cast<std::size_t>(values.size());
  BollingerBands bands{math::unobserved(rows), math::unobserved(rows), math::unobserved(rows)};
  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_) {
      continue;
    }
    bands.middle_[row] = spread.mean_;
    bands.upper_[row] = spread.mean_ + multiplier * spread.deviation_;
    bands.lower_[row] = spread.mean_ - multiplier * spread.deviation_;
  }
  return bands;
}

inline BollingerBands bollinger(
  math::SeriesView values,
  std::size_t window = 20,
  Scalar multiplier = Scalar{2}
) {
  return bollinger(values, window, multiplier, window);
}

inline math::Vector rsi(math::SeriesView prices, std::size_t period = 14) {
  if (period == 0) {
    throw std::invalid_argument("rsi period must be at least one row");
  }
  const auto rows = static_cast<std::size_t>(prices.size());
  math::Vector gains = math::unobserved(rows);
  math::Vector losses = math::unobserved(rows);
  for (math::Index row = 1; row < prices.size(); ++row) {
    if (!observed(prices[row - 1]) || !observed(prices[row])) {
      continue;
    }
    const Scalar change = prices[row] - prices[row - 1];
    gains[row] = std::max(change, Scalar{0});
    losses[row] = std::max(-change, Scalar{0});
  }
 
  const math::Vector average_gain = indicator_detail::wilder(gains, period);
  const math::Vector average_loss = indicator_detail::wilder(losses, period);
 
  math::Vector result = math::unobserved(rows);
  for (math::Index row = 0; row < prices.size(); ++row) {
    if (!observed(average_gain[row]) || !observed(average_loss[row])) {
      continue;
    }
    if (average_loss[row] == Scalar{0}) {
      result[row] = average_gain[row] == Scalar{0} ? Scalar{50} : Scalar{100};
      continue;
    }
    const Scalar strength = average_gain[row] / average_loss[row];
    result[row] = Scalar{100} - Scalar{100} / (Scalar{1} + strength);
  }
  return result;
}

inline math::Vector true_range(
  math::SeriesView highs,
  math::SeriesView lows,
  math::SeriesView closes
) {
  indicator_detail::require_same_rows(highs, lows);
  indicator_detail::require_same_rows(highs, closes);
 
  math::Vector result = math::unobserved(static_cast<std::size_t>(highs.size()));
  for (math::Index row = 0; row < highs.size(); ++row) {
    if (!observed(highs[row]) || !observed(lows[row])) {
      continue;
    }
    Scalar range = highs[row] - lows[row];
    if (row > 0 && observed(closes[row - 1])) {
      const Scalar previous = closes[row - 1];
      range = std::max(range, std::abs(highs[row] - previous));
      range = std::max(range, std::abs(lows[row] - previous));
    }
    result[row] = range;
  }
  return result;
}

inline math::Vector atr(
  math::SeriesView highs,
  math::SeriesView lows,
  math::SeriesView closes,
  std::size_t period = 14
) {
  if (period == 0) {
    throw std::invalid_argument("atr period must be at least one bar");
  }
  return indicator_detail::wilder(true_range(highs, lows, closes), period);
}

struct MacdLines {
  math::Vector macd_;
  math::Vector signal_;
  math::Vector histogram_;
};

inline MacdLines macd(
  math::SeriesView values,
  std::size_t fast = 12,
  std::size_t slow = 26,
  std::size_t signal_span = 9
) {
  if (fast >= slow) {
    throw std::invalid_argument("macd fast span must be shorter than the slow span");
  }
  const math::Vector quick = ema(values, alpha_from_span(fast));
  const math::Vector slower = ema(values, alpha_from_span(slow));
 
  const auto rows = static_cast<std::size_t>(values.size());
  MacdLines lines{math::unobserved(rows), math::unobserved(rows), math::unobserved(rows)};
  for (math::Index row = 0; row < values.size(); ++row) {
    if (observed(quick[row]) && observed(slower[row])) {
      lines.macd_[row] = quick[row] - slower[row];
    }
  }
  lines.signal_ = ema(lines.macd_, alpha_from_span(signal_span));
  for (math::Index row = 0; row < values.size(); ++row) {
    if (observed(lines.macd_[row]) && observed(lines.signal_[row])) {
      lines.histogram_[row] = lines.macd_[row] - lines.signal_[row];
    }
  }
  return lines;
}

inline math::Vector rolling_cov(
  math::SeriesView left,
  math::SeriesView right,
  std::size_t window,
  std::size_t min_observed
) {
  indicator_detail::require_same_rows(left, right);
  rolling_detail::check(window, min_observed);
  const auto paired = indicator_detail::pair_series(left, right);
  math::Vector result = math::unobserved(static_cast<std::size_t>(left.size()));
  for (math::Index row = 0; row < left.size(); ++row) {
    const math::Index first = rolling_detail::window_start(row, window);
    const math::Index length = row - first + 1;
    const auto reduced = indicator_detail::reduce_pair(paired, first, length, min_observed);
    if (!reduced.reported_) {
      continue;
    }
    const math::Vector left_centred = indicator_detail::centred(
      paired.left_, paired.present_, first, length, reduced.mean_left_
    );
    const math::Vector right_centred = indicator_detail::centred(
      paired.right_, paired.present_, first, length, reduced.mean_right_
    );
    result[row] = left_centred.dot(right_centred) / static_cast<Scalar>(reduced.count_ - 1);
  }
  return result;
}

inline math::Vector rolling_cov(
  math::SeriesView left,
  math::SeriesView right,
  std::size_t window
) {
  return rolling_cov(left, right, window, window);
}

inline math::Vector rolling_corr(
  math::SeriesView left,
  math::SeriesView right,
  std::size_t window,
  std::size_t min_observed
) {
  indicator_detail::require_same_rows(left, right);
  rolling_detail::check(window, min_observed);
  const auto paired = indicator_detail::pair_series(left, right);
  math::Vector result = math::unobserved(static_cast<std::size_t>(left.size()));
  for (math::Index row = 0; row < left.size(); ++row) {
    const math::Index first = rolling_detail::window_start(row, window);
    const math::Index length = row - first + 1;
    const auto reduced = indicator_detail::reduce_pair(paired, first, length, min_observed);
    if (!reduced.reported_) {
      continue;
    }
    const math::Vector left_centred = indicator_detail::centred(
      paired.left_, paired.present_, first, length, reduced.mean_left_
    );
    const math::Vector right_centred = indicator_detail::centred(
      paired.right_, paired.present_, first, length, reduced.mean_right_
    );
    const Scalar squared_left = left_centred.squaredNorm();
    const Scalar squared_right = right_centred.squaredNorm();
    const auto spread = [&](Scalar squared) {
      return std::sqrt(squared / static_cast<Scalar>(reduced.count_ - 1));
    };
    if (squared_left == Scalar{0} || squared_right == Scalar{0} ||
        constant_detail::rounded_constant(
          spread(squared_left), reduced.mean_left_, reduced.count_, paired.left_, paired.present_, first, length
        ) ||
        constant_detail::rounded_constant(
          spread(squared_right), reduced.mean_right_, reduced.count_, paired.right_, paired.present_, first, length
        )) {
      // A constant series. Correlation with it is undefined, not zero.
      result[row] = math::undefined();
      continue;
    }
    // The root of the product, which keeps a collinear pair at exactly 1 or -1; two roots only where the
    // product of two finite squared sums underflows (tiny dispersions) or overflows (large ones).
    const Scalar product = squared_left * squared_right;
    const Scalar denominator =
      std::isnormal(product) ? std::sqrt(product) : std::sqrt(squared_left) * std::sqrt(squared_right);
    result[row] = left_centred.dot(right_centred) / denominator;
  }
  return result;
}

inline math::Vector rolling_corr(
  math::SeriesView left,
  math::SeriesView right,
  std::size_t window
) {
  return rolling_corr(left, right, window, window);
}

inline math::Vector beta(
  math::SeriesView asset,
  math::SeriesView market,
  std::size_t window,
  std::size_t min_observed
) {
  indicator_detail::require_same_rows(asset, market);
  rolling_detail::check(window, min_observed);
  const auto paired = indicator_detail::pair_series(asset, market);
  math::Vector result = math::unobserved(static_cast<std::size_t>(asset.size()));
  for (math::Index row = 0; row < asset.size(); ++row) {
    const math::Index first = rolling_detail::window_start(row, window);
    const math::Index length = row - first + 1;
    const auto reduced = indicator_detail::reduce_pair(paired, first, length, min_observed);
    if (!reduced.reported_) {
      continue;
    }
    const math::Vector market_centred = indicator_detail::centred(
      paired.right_, paired.present_, first, length, reduced.mean_right_
    );
    const Scalar squared_market = market_centred.squaredNorm();
    if (squared_market == Scalar{0} ||
        constant_detail::rounded_constant(
          std::sqrt(squared_market / static_cast<Scalar>(reduced.count_ - 1)),
          reduced.mean_right_,
          reduced.count_,
          paired.right_,
          paired.present_,
          first,
          length
        )) {
      // No market movement to attribute the asset's movement to.
      result[row] = math::undefined();
      continue;
    }
    const math::Vector asset_centred = indicator_detail::centred(
      paired.left_, paired.present_, first, length, reduced.mean_left_
    );
    result[row] = asset_centred.dot(market_centred) / squared_market;
  }
  return result;
}

inline math::Vector beta(
  math::SeriesView asset,
  math::SeriesView market,
  std::size_t window
) {
  return beta(asset, market, window, window);
}
 
} // namespace analytics
} // namespace algotradelib