// RAID vs R (highfrequency 1.0.3, rugarch 1.5.6) on the committed fixtures — values from scripts/reference_r.R.
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <nlohmann/json.hpp>
#include <fstream>
#include "raid/data/io.hpp"
#include "raid/data/returns.hpp"
#include "raid/estimators/realized.hpp"
#include "raid/estimators/har.hpp"
#include "raid/models/garch.hpp"
using namespace raid; using json = nlohmann::json; using Catch::Matchers::WithinRel; using Catch::Matchers::WithinAbs;

static json load() { std::ifstream f(std::string(RAID_SOURCE_DIR) + "/datasets/fixtures/reference_r.json"); REQUIRE(f.good()); json j; f >> j; return j; }

TEST_CASE("Realized measures match highfrequency (rCov, rBPCov, rSVar, rQuar)", "[reference][r][realized]") {
  auto R = load(); auto bars = data::read_bars_csv(std::string(RAID_SOURCE_DIR) + "/datasets/fixtures/synthetic_bars_1m.csv"); MarketCalendar cal;
  for (auto& e : R["realized"]) {
    const std::string sym = e["symbol"], sess = e["session"]; auto b = bars.at(sym); b.normalize(); auto reg = b.filter_session(cal, SessionPolicy::regular); auto ir = data::intraday_returns(reg);
    usize k = 0; while (k < ir.session_dates.size() && ir.session_dates[k].iso() != sess) ++k; REQUIRE(k < ir.session_dates.size());
    est::RealizedOptions o; o.plan.backend = Backend::cpu; o.bpv_finite_sample = false; auto m = est::realized_measures(ir, o);   // highfrequency's rBPCov with makeReturns=FALSE applies no M/(M−1) factor
    INFO(sym << " " << sess); CHECK(m.n_returns[k] == e["n_returns"].get<usize>());
    CHECK_THAT(m.rv[k], WithinRel(e["rv"].get<double>(), 1e-10)); CHECK_THAT(m.bpv[k], WithinRel(e["bpv"].get<double>(), 1e-10));
    CHECK_THAT(m.rsv_neg[k], WithinRel(e["rsv"][0].get<double>(), 1e-10)); CHECK_THAT(m.rsv_pos[k], WithinRel(e["rsv"][1].get<double>(), 1e-10));
    // highfrequency::rQuar uses N = (number of returns + 1)/3 (number of prices); RAID uses M/3 (BPQ 2016, BNS 2002): convert
    const double M = static_cast<double>(m.n_returns[k]); CHECK_THAT(m.rq[k] * (M + 1.0) / M, WithinRel(e["rq"].get<double>(), 1e-10));
  }
}
TEST_CASE("HAR-RV coefficients match highfrequency::HARmodel (h = 1)", "[reference][r][har]") {
  auto R = load(); std::vector<double> rv, rq; { std::ifstream f(std::string(RAID_SOURCE_DIR) + "/datasets/fixtures/har_rv.csv"); std::string line; std::getline(f, line); while (std::getline(f, line)) { auto c = line.find(','); rv.push_back(std::stod(line.substr(0, c))); rq.push_back(std::stod(line.substr(c + 1))); } }
  est::HarSpec s; s.horizon = 1; auto e = est::har_fit(rv, s); auto h1 = R["har_h1"];
  CHECK(e.n_obs == h1["n"].get<usize>()); for (int j = 0; j < 4; ++j) CHECK_THAT(e.coefficients[j], WithinRel(h1["coef"][j].get<double>(), 1e-8)); CHECK_THAT(e.loss, WithinRel(h1["ssr"].get<double>(), 1e-8));
  // h = 5: highfrequency uses the h-day *average* RV as the dependent variable (same as our HarTarget::average)
  s.horizon = 5; auto e5 = est::har_fit(rv, s); auto h5 = R["har_h5"]; INFO("n ours " << e5.n_obs << " R " << h5["n"].get<usize>());
  if (e5.n_obs == h5["n"].get<usize>()) { for (int j = 0; j < 4; ++j) CHECK_THAT(e5.coefficients[j], WithinRel(h5["coef"][j].get<double>(), 1e-6)); }
  else WARN("HARmodel h=5 sample size differs (" + std::to_string(h5["n"].get<usize>()) + " vs " + std::to_string(e5.n_obs) + "): alignment convention documented in the registry");
  // HARQ: highfrequency's RQ1 term is sqrt(RQ_t) * RV_t (BPQ 2016); coefficient scale must agree
  est::HarSpec q; q.harq = true; auto eq = est::har_fit(rv, q, rq); auto hq = R["harq_h1"]; INFO("harq n ours " << eq.n_obs << " R " << hq["n"].get<usize>());
  if (eq.n_obs == hq["n"].get<usize>()) for (int j = 0; j < 5; ++j) CHECK_THAT(eq.coefficients[j], WithinRel(hq["coef"][j].get<double>(), 1e-6));
}
TEST_CASE("GARCH family matches rugarch within tolerance (parameter mapping documented)", "[reference][r][garch]") {
  auto R = load(); std::vector<double> r; { std::ifstream f(std::string(RAID_SOURCE_DIR) + "/datasets/fixtures/garch_returns.csv"); std::string line; std::getline(f, line); while (std::getline(f, line)) r.push_back(std::stod(line)); }
  using namespace models;
  auto check = [&](const char* key, GarchModel m, GarchDist d, std::vector<std::pair<int, const char*>> map) {
    GarchFitOptions o; o.model = m; o.dist = d; auto fit = garch_fit(r, o); INFO(key); REQUIRE(fit.ok()); auto g = R["garch"][key];
    // rugarch initialises σ² with the sample variance (as RAID) → log-likelihoods should agree to 1e-4 relative
    CHECK_THAT(-fit.loss, WithinRel(g["loglik"].get<double>(), 2e-4));
    for (auto& [j, name] : map) { const double ours = fit.coefficients[static_cast<usize>(j)], theirs = g["coef"][name].get<double>(); INFO(name << " ours " << ours << " rugarch " << theirs); CHECK((std::abs(ours - theirs) <= 0.05 * std::abs(theirs) + 0.01)); }
    return fit; };
  check("sGARCH.norm", GarchModel::garch, GarchDist::gaussian, {{0, "mu"}, {1, "omega"}, {2, "alpha1"}, {3, "beta1"}});
  check("gjrGARCH.norm", GarchModel::gjr, GarchDist::gaussian, {{0, "mu"}, {1, "omega"}, {2, "alpha1"}, {3, "gamma1"}, {4, "beta1"}});
  // rugarch eGARCH: log σ² = ω + α z + γ(|z| − E|z|) + β log σ²  ⇒ rugarch alpha1 = RAID gamma (sign), rugarch gamma1 = RAID alpha (magnitude)
  check("eGARCH.norm", GarchModel::egarch, GarchDist::gaussian, {{0, "mu"}, {1, "omega"}, {2, "gamma1"}, {3, "alpha1"}, {4, "beta1"}});
  auto t = check("sGARCH.std", GarchModel::garch, GarchDist::student_t, {{0, "mu"}, {1, "omega"}, {2, "alpha1"}, {3, "beta1"}});
  CHECK(t.coefficients.back() > 30.0);   // rugarch shape → 100 (its upper bound); RAID ν → its own upper bound (500): Gaussian data, ν unidentified
}