// Reference tests: RAID vs established software on the committed fixtures (values produced by scripts/make_fixtures.py
// with statsmodels 0.15 / scipy 1.18 / arch 8.0 / exchange_calendars 4.13; versions recorded in reference_values.json).
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <nlohmann/json.hpp>
#include <fstream>
#include <set>
#include "raid/estimators/ols.hpp"
#include "raid/estimators/realized.hpp"
#include "raid/models/garch.hpp"
#include "raid/data/calendar.hpp"
#include "raid/data/io.hpp"
#include "raid/data/returns.hpp"
using namespace raid; using json = nlohmann::json; using Catch::Matchers::WithinRel; using Catch::Matchers::WithinAbs;

static json load_ref() { std::ifstream f(std::string(RAID_SOURCE_DIR) + "/datasets/fixtures/reference_values.json"); REQUIRE(f.good()); json j; f >> j; return j; }
static void load_csv(const std::string& p, std::vector<std::vector<double>>& rows) { std::ifstream f(p); REQUIRE(f.good()); std::string line; std::getline(f, line); while (std::getline(f, line)) { if (line.empty()) continue; std::vector<double> r; usize s = 0; for (usize i = 0; i <= line.size(); ++i) if (i == line.size() || line[i] == ',') { r.push_back(std::stod(line.substr(s, i - s))); s = i + 1; } rows.push_back(r); } }

TEST_CASE("OLS and robust standard errors match statsmodels", "[reference][ols]") {
  auto ref = load_ref()["ols"]; std::vector<std::vector<double>> rows; load_csv(std::string(RAID_SOURCE_DIR) + "/datasets/fixtures/ols_design.csv", rows);
  const Eigen::Index n = static_cast<Eigen::Index>(rows.size()); Eigen::MatrixXd X(n, 3); Eigen::VectorXd y(n);
  for (Eigen::Index i = 0; i < n; ++i) { y[i] = rows[i][0]; X(i, 0) = 1; X(i, 1) = rows[i][1]; X(i, 2) = rows[i][2]; }
  auto r = est::ols(X, y);
  for (int j = 0; j < 3; ++j) { CHECK_THAT(r.estimate.coefficients[j], WithinRel(ref["params"][j].get<double>(), 1e-10)); CHECK_THAT(r.estimate.inference->std_errors[j], WithinRel(ref["bse"][j].get<double>(), 1e-9)); CHECK_THAT(r.estimate.inference->p_values[j], WithinAbs(ref["pvalues"][j].get<double>(), 1e-9)); }
  CHECK_THAT(r.r2, WithinRel(ref["rsquared"].get<double>(), 1e-10)); CHECK_THAT(r.estimate.loss, WithinRel(ref["ssr"].get<double>(), 1e-10));
  est::OlsOptions o; o.cov.type = stats::CovType::hc0; auto h0 = est::ols(X, y, o); o.cov.type = stats::CovType::hc3; auto h3 = est::ols(X, y, o);
  o.cov.type = stats::CovType::hac_bartlett; o.cov.hac_lags = 5; auto hac = est::ols(X, y, o);
  std::vector<i64> g(static_cast<usize>(n)); for (Eigen::Index i = 0; i < n; ++i) g[static_cast<usize>(i)] = i / 25; o.cov.type = stats::CovType::cluster; o.cov.clusters = g; auto cl = est::ols(X, y, o);
  for (int j = 0; j < 3; ++j) { CHECK_THAT(h0.estimate.inference->std_errors[j], WithinRel(ref["hc0_bse"][j].get<double>(), 1e-9)); CHECK_THAT(h3.estimate.inference->std_errors[j], WithinRel(ref["hc3_bse"][j].get<double>(), 1e-9)); CHECK_THAT(hac.estimate.inference->std_errors[j], WithinRel(ref["hac5_bse"][j].get<double>(), 1e-9)); CHECK_THAT(cl.estimate.inference->std_errors[j], WithinRel(ref["cluster_bse"][j].get<double>(), 1e-9)); }
  CHECK_THAT(est::t_pvalue_two_sided(1.96, 0), WithinAbs(ref["t_pvalues"]["1.96_inf"].get<double>(), 1e-12));
  CHECK_THAT(est::t_pvalue_two_sided(2.0, 10), WithinAbs(ref["t_pvalues"]["2.0_10"].get<double>(), 1e-12));
}

TEST_CASE("XNYS calendar matches exchange_calendars", "[reference][calendar]") {
  auto ref = load_ref()["xnys"]; if (ref.contains("error")) { WARN("exchange_calendars unavailable when fixtures were generated: " + ref["error"].get<std::string>()); return; }
  MarketCalendar cal; std::set<std::string> ours; for (auto d : cal.sessions(Date::from_ymd(2000, 1, 1), Date::from_ymd(2026, 12, 31))) ours.insert(d.iso());
  CHECK(ours.size() == ref["n_sessions_2000_2026"].get<usize>());
  std::set<std::string> theirs2024; for (auto& s : ref["sessions_2024"]) theirs2024.insert(s.get<std::string>());
  std::set<std::string> ours2024; for (auto d : cal.sessions(Date::from_ymd(2024, 1, 1), Date::from_ymd(2024, 12, 31))) ours2024.insert(d.iso());
  CHECK(ours2024 == theirs2024);
  usize mismatches = 0; for (auto& h : ref["holidays_2000_2026"]) if (ours.count(h.get<std::string>())) { ++mismatches; UNSCOPED_INFO("RAID trades on " << h.get<std::string>() << " but exchange_calendars does not"); }
  CHECK(mismatches == 0);
  if (!ref["early_closes_2000_2026"].is_null()) { usize bad = 0; for (auto& e : ref["early_closes_2000_2026"]) { Date d = date_from_iso(e.get<std::string>()); if (!cal.is_early_close(d)) { ++bad; UNSCOPED_INFO("missing early close " << e.get<std::string>()); } } CHECK(bad == 0);
    usize extra = 0; for (auto d : cal.sessions(Date::from_ymd(2000, 1, 1), Date::from_ymd(2026, 12, 31))) if (cal.is_early_close(d)) { bool found = false; for (auto& e : ref["early_closes_2000_2026"]) if (e.get<std::string>() == d.iso()) found = true; if (!found) { ++extra; UNSCOPED_INFO("extra early close " << d.iso()); } } CHECK(extra == 0); }
}

TEST_CASE("Realized measures match the fixture truth (numpy)", "[reference][realized]") {
  std::ifstream f(std::string(RAID_SOURCE_DIR) + "/datasets/fixtures/synthetic_bars_1m.expected.json"); json truth; f >> truth;
  auto bars = data::read_bars_csv(std::string(RAID_SOURCE_DIR) + "/datasets/fixtures/synthetic_bars_1m.csv"); MarketCalendar cal;
  for (auto& [sym, b0] : bars) { auto b = b0; b.normalize(); auto reg = b.filter_session(cal, SessionPolicy::regular); auto ir = data::intraday_returns(reg); est::RealizedOptions o; o.plan.backend = Backend::cpu; auto m = est::realized_measures(ir, o);
    REQUIRE(m.sessions.size() == 5);
    for (usize k = 0; k < 5; ++k) { auto t = truth[sym][m.sessions[k].iso()]; CHECK(m.n_returns[k] + 1 == t["n_regular_bars"].get<usize>()); CHECK(cal.expected_regular_bars(m.sessions[k], Interval::m1) == t["expected_bars"].get<usize>());
      CHECK_THAT(m.rv[k], WithinRel(t["rv"].get<double>(), 1e-12)); CHECK_THAT(m.bpv[k], WithinRel(t["bpv"].get<double>(), 1e-12)); CHECK_THAT(m.rsv_neg[k], WithinRel(t["rsv_neg"].get<double>(), 1e-12)); CHECK_THAT(m.rsv_pos[k], WithinRel(t["rsv_pos"].get<double>(), 1e-12)); CHECK_THAT(m.rq[k], WithinRel(t["rq"].get<double>(), 1e-12)); } }
}

TEST_CASE("GARCH-family QML matches the arch package within tolerance", "[reference][garch]") {
  auto ref = load_ref()["garch"]; if (ref.contains("error")) { WARN("arch unavailable: " + ref["error"].get<std::string>()); return; }
  std::vector<std::vector<double>> rows; load_csv(std::string(RAID_SOURCE_DIR) + "/datasets/fixtures/garch_returns.csv", rows); std::vector<double> r; for (auto& x : rows) r.push_back(x[0]);
  struct Case { const char* key; models::GarchModel m; models::GarchDist d; std::vector<const char*> arch_names; };
  const std::vector<Case> cases = {{"garch.normal", models::GarchModel::garch, models::GarchDist::gaussian, {"mu", "omega", "alpha[1]", "beta[1]"}},
                                   {"gjr.normal", models::GarchModel::gjr, models::GarchDist::gaussian, {"mu", "omega", "alpha[1]", "gamma[1]", "beta[1]"}},
                                   {"garch.t", models::GarchModel::garch, models::GarchDist::student_t, {"mu", "omega", "alpha[1]", "beta[1]", "nu"}},
                                   {"egarch.normal", models::GarchModel::egarch, models::GarchDist::gaussian, {"mu", "omega", "alpha[1]", "gamma[1]", "beta[1]"}}};
  for (auto& c : cases) {
    auto a = ref[c.key]; models::GarchFitOptions o; o.model = c.m; o.dist = c.d; o.forecast_horizon = 5; auto fit = models::garch_fit(r, o);
    INFO(c.key); REQUIRE(fit.ok());
    // Both maximise the same likelihood up to the initial-variance convention: log-likelihoods within 0.05 % and
    // parameters within 5 % relative or 0.01 absolute (documented tolerance, see docs/methods/registry.md).
    CHECK_THAT(-fit.loss, WithinRel(a["loglik"].get<double>(), 5e-4));
    for (usize j = 0; j < c.arch_names.size(); ++j) { const double ours = fit.coefficients[j], theirs = a["params"][c.arch_names[j]].get<double>();
      if (std::string(c.arch_names[j]) == "nu" && theirs > 100) continue;  // ν unidentified for Gaussian data
      INFO(c.arch_names[j] << " ours=" << ours << " arch=" << theirs); CHECK((std::abs(ours - theirs) <= 0.05 * std::abs(theirs) + 0.01)); }
    if (c.m != models::GarchModel::egarch) { auto fc = fit.extra["forecast_variance"].get<std::vector<double>>(); for (int h = 0; h < 5; ++h) CHECK_THAT(fc[static_cast<usize>(h)], WithinRel(a["forecast_h5"][h].get<double>(), 0.05)); }
  }
}