1// Reference tests: RAID vs established software on the committed fixtures (values produced by scripts/make_fixtures.py2// with statsmodels 0.15 / scipy 1.18 / arch 8.0 / exchange_calendars 4.13; versions recorded in reference_values.json).3#include <catch2/catch_test_macros.hpp>4#include <catch2/matchers/catch_matchers_floating_point.hpp>5#include <nlohmann/json.hpp>6#include <fstream>7#include <set>8#include "raid/estimators/ols.hpp"9#include "raid/estimators/realized.hpp"10#include "raid/models/garch.hpp"11#include "raid/data/calendar.hpp"12#include "raid/data/io.hpp"13#include "raid/data/returns.hpp"14using namespace raid; using json = nlohmann::json; using Catch::Matchers::WithinRel; using Catch::Matchers::WithinAbs;1516static 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; }17static 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); } }1819TEST_CASE("OLS and robust standard errors match statsmodels", "[reference][ols]") {20 auto ref = load_ref()["ols"]; std::vector<std::vector<double>> rows; load_csv(std::string(RAID_SOURCE_DIR) + "/datasets/fixtures/ols_design.csv", rows);21 const Eigen::Index n = static_cast<Eigen::Index>(rows.size()); Eigen::MatrixXd X(n, 3); Eigen::VectorXd y(n);22 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]; }23 auto r = est::ols(X, y);24 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)); }25 CHECK_THAT(r.r2, WithinRel(ref["rsquared"].get<double>(), 1e-10)); CHECK_THAT(r.estimate.loss, WithinRel(ref["ssr"].get<double>(), 1e-10));26 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);27 o.cov.type = stats::CovType::hac_bartlett; o.cov.hac_lags = 5; auto hac = est::ols(X, y, o);28 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);29 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)); }30 CHECK_THAT(est::t_pvalue_two_sided(1.96, 0), WithinAbs(ref["t_pvalues"]["1.96_inf"].get<double>(), 1e-12));31 CHECK_THAT(est::t_pvalue_two_sided(2.0, 10), WithinAbs(ref["t_pvalues"]["2.0_10"].get<double>(), 1e-12));32}3334TEST_CASE("XNYS calendar matches exchange_calendars", "[reference][calendar]") {35 auto ref = load_ref()["xnys"]; if (ref.contains("error")) { WARN("exchange_calendars unavailable when fixtures were generated: " + ref["error"].get<std::string>()); return; }36 MarketCalendar cal; std::set<std::string> ours; for (auto d : cal.sessions(Date::from_ymd(2000, 1, 3), Date::from_ymd(2026, 12, 31))) ours.insert(d.iso());37 CHECK(ours.size() == ref["n_sessions_2000_2026"].get<usize>());38 std::set<std::string> theirs2024; for (auto& s : ref["sessions_2024"]) theirs2024.insert(s.get<std::string>());39 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());40 CHECK(ours2024 == theirs2024);41 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"); }42 CHECK(mismatches == 0);43 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);44 usize extra = 0; for (auto d : cal.sessions(Date::from_ymd(2000, 1, 3), 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); }45}4647TEST_CASE("Realized measures match the fixture truth (numpy)", "[reference][realized]") {48 std::ifstream f(std::string(RAID_SOURCE_DIR) + "/datasets/fixtures/synthetic_bars_1m.expected.json"); json truth; f >> truth;49 auto bars = data::read_bars_csv(std::string(RAID_SOURCE_DIR) + "/datasets/fixtures/synthetic_bars_1m.csv"); MarketCalendar cal;50 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);51 REQUIRE(m.sessions.size() == 5);52 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>());53 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)); } }54}5556TEST_CASE("GARCH-family QML matches the arch package within tolerance", "[reference][garch]") {57 auto ref = load_ref()["garch"]; if (ref.contains("error")) { WARN("arch unavailable: " + ref["error"].get<std::string>()); return; }58 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]);59 struct Case { const char* key; models::GarchModel m; models::GarchDist d; std::vector<const char*> arch_names; };60 const std::vector<Case> cases = {{"garch.normal", models::GarchModel::garch, models::GarchDist::gaussian, {"mu", "omega", "alpha[1]", "beta[1]"}},61 {"gjr.normal", models::GarchModel::gjr, models::GarchDist::gaussian, {"mu", "omega", "alpha[1]", "gamma[1]", "beta[1]"}},62 {"garch.t", models::GarchModel::garch, models::GarchDist::student_t, {"mu", "omega", "alpha[1]", "beta[1]", "nu"}},63 {"egarch.normal", models::GarchModel::egarch, models::GarchDist::gaussian, {"mu", "omega", "alpha[1]", "gamma[1]", "beta[1]"}}};64 for (auto& c : cases) {65 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);66 INFO(c.key); REQUIRE(fit.ok());67 // Both maximise the same likelihood up to the initial-variance convention: log-likelihoods within 0.05 % and68 // parameters within 5 % relative or 0.01 absolute (documented tolerance, see docs/methods/registry.md).69 CHECK_THAT(-fit.loss, WithinRel(a["loglik"].get<double>(), 5e-4));70 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>();71 if (std::string(c.arch_names[j]) == "nu" && theirs > 100) continue; // ν unidentified for Gaussian data72 INFO(c.arch_names[j] << " ours=" << ours << " arch=" << theirs); CHECK((std::abs(ours - theirs) <= 0.05 * std::abs(theirs) + 0.01)); }73 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)); }74 }75}76