metrics.hpp

Go to the documentation of this file.

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

 
namespace metric_detail {

struct Summary {
  std::size_t count_{0};
  Scalar mean_{0};
  Scalar squared_{0};
  // Whether at least two observations are present and all equal, as far as rounding can hide it
  // (`constant_detail::rounded_constant`); an exactly zero deviation is the caller's own test.
  bool constant_{false};
};

inline Summary summarize(math::SeriesView values) {
  const auto parts = math::observed_parts(values);
  Summary result;
  const Scalar count = parts.present_.sum();
  result.count_ = static_cast<std::size_t>(count);
  if (result.count_ == 0) {
    return result;
  }
  result.mean_ = parts.values_.sum() / count;
  result.squared_ = (parts.values_ - result.mean_ * parts.present_).squaredNorm();
  if (result.count_ >= 2) {
    const Scalar deviation = std::sqrt(result.squared_ / static_cast<Scalar>(result.count_ - 1));
    result.constant_ = constant_detail::rounded_constant(
      deviation, result.mean_, result.count_, parts.values_, parts.present_, 0, parts.values_.size()
    );
  }
  return result;
}
 
} // namespace metric_detail

struct MaxDrawdown {
  Scalar magnitude_{0};
  math::Index peak_row_{-1};
  math::Index trough_row_{-1};
};

inline MaxDrawdown max_drawdown(math::SeriesView levels) {
  MaxDrawdown result;
  result.magnitude_ = missing();
  Scalar peak{0};
  math::Index peak_row = -1;
  std::size_t observations = 0;
  for (math::Index row = 0; row < levels.size(); ++row) {
    if (!observed(levels[row])) {
      continue;
    }
    ++observations;
    if (peak_row < 0 || peak < levels[row]) {
      peak = levels[row];
      peak_row = row;
    }
    if (peak <= Scalar{0}) {
      continue;
    }
    const Scalar decline = (peak - levels[row]) / peak;
    if (!observed(result.magnitude_) || result.magnitude_ < decline) {
      result.magnitude_ = decline;
      result.peak_row_ = peak_row;
      result.trough_row_ = row;
    }
  }
  // Levels were present but no ratio could be formed from them -- every running peak was at or
  // below zero. That is "undefined", not "no data".
  if (observations > 0 && !observed(result.magnitude_)) {
    result.magnitude_ = math::undefined();
  }
  return result;
}

inline Scalar sharpe(math::SeriesView period_returns, Scalar risk_free = Scalar{0}) {
  const auto summary = metric_detail::summarize(period_returns);
  if (summary.count_ < 2) {
    return missing();
  }
  const Scalar deviation =
    std::sqrt(summary.squared_ / static_cast<Scalar>(summary.count_ - 1));
  if (deviation == Scalar{0} || summary.constant_) {
    // Returns were present and never varied. `undefined()` rather than `missing()`, because a
    // caller ranking strategies should be able to tell this from "not enough history".
    return math::undefined();
  }
  return (summary.mean_ - risk_free) / deviation;
}

inline Scalar sortino(math::SeriesView period_returns, Scalar target = Scalar{0}) {
  std::size_t count = 0;
  std::size_t below = 0;
  Scalar sum{0};
  Scalar shortfall{0};
  for (math::Index row = 0; row < period_returns.size(); ++row) {
    const Scalar value = period_returns[row];
    if (!observed(value)) {
      continue;
    }
    ++count;
    sum += value;
    if (value < target) {
      ++below;
      const Scalar gap = target - value;
      shortfall += gap * gap;
    }
  }
  if (count < 2) {
    return missing();
  }
  if (below == 0) {
    // Nothing fell short of the target. There is no downside risk to divide by, and an infinity
    // would survive every later comparison as the best strategy ever measured. `undefined()`
    // says "computed, and does not exist" rather than "no data", so a caller ranking strategies
    // can treat the two differently.
    return math::undefined();
  }
  // Divided by the whole count, not by `below`. That is the definition, and the difference
  // is a systematic flattery of strategies with rare losses.
  const Scalar downside = std::sqrt(shortfall / static_cast<Scalar>(count));
  if (downside == Scalar{0}) {
    return math::undefined();
  }
  return (sum / static_cast<Scalar>(count) - target) / downside;
}

inline Scalar annualize_ratio(Scalar ratio, Scalar periods_per_year) {
  if (!(periods_per_year > Scalar{0})) {
    throw std::invalid_argument("periods per year must be positive");
  }
  return observed(ratio) ? ratio * std::sqrt(periods_per_year) : ratio;
}

enum class TurnoverStart {
  FROM_CASH,
  SKIP,
};

inline math::Vector one_way_turnover(
  math::MatrixView weights,
  TurnoverStart start = TurnoverStart::FROM_CASH
) {
  // Masked once for the whole panel, so each period is a plain vectorized row difference.
  const math::Matrix held =
    (weights.array() == weights.array()).select(weights, Scalar{0});
  math::Vector result = math::unobserved(static_cast<std::size_t>(weights.rows()));
  for (math::Index row = 0; row < held.rows(); ++row) {
    if (row == 0) {
      if (start == TurnoverStart::SKIP) {
        continue;
      }
      result[row] = Scalar{0.5} * held.row(0).cwiseAbs().sum();
      continue;
    }
    result[row] = Scalar{0.5} * (held.row(row) - held.row(row - 1)).cwiseAbs().sum();
  }
  return result;
}

inline math::Vector two_way_turnover(
  math::MatrixView weights,
  TurnoverStart start = TurnoverStart::FROM_CASH
) {
  return Scalar{2} * one_way_turnover(weights, start);
}
 
} // namespace analytics
} // namespace algotradelib