1#include <catch2/catch_test_macros.hpp>2#include <catch2/matchers/catch_matchers_floating_point.hpp>3#include <nlohmann/json.hpp>4#include <fstream>5#include "raid/models/var.hpp"6#include "raid/core/rng.hpp"7using namespace raid; using json = nlohmann::json; using Catch::Matchers::WithinRel; using Catch::Matchers::WithinAbs;89static Eigen::MatrixXd load_var_csv(int& N) { std::ifstream f(std::string(RAID_SOURCE_DIR) + "/datasets/fixtures/var_panel.csv"); REQUIRE(f.good()); std::string line; std::getline(f, line); std::vector<std::vector<double>> rows; while (std::getline(f, line)) { 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); } N = static_cast<int>(rows[0].size()); Eigen::MatrixXd Y(rows.size(), N); for (usize i = 0; i < rows.size(); ++i) for (int j = 0; j < N; ++j) Y(static_cast<Eigen::Index>(i), j) = rows[i][static_cast<usize>(j)]; return Y; }1011TEST_CASE("VAR(2) coefficients, Σ, FEVD and generalized FEVD match statsmodels", "[reference][var]") {12 std::ifstream f(std::string(RAID_SOURCE_DIR) + "/datasets/fixtures/reference_values.json"); json ref; f >> ref; auto v = ref["var"]; int N; auto Y = load_var_csv(N);13 auto fit = models::var_fit(Y, 2); CHECK(fit.n_obs == v["nobs"].get<usize>()); CHECK(fit.stable);14 // statsmodels params: (1 + N p) × N, rows = [const, L1.y1..L1.yN, L2...], our coef is N × (1+Np)15 for (int i = 0; i < N; ++i) for (int c = 0; c < 1 + 2 * N; ++c) CHECK_THAT(fit.coef(i, c), WithinAbs(v["params"][c][i].get<double>(), 1e-9));16 for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) CHECK_THAT(fit.sigma(i, j), WithinRel(v["sigma_u"][i][j].get<double>(), 1e-9));17 CHECK_THAT(fit.loglik, WithinRel(v["llf"].get<double>(), 1e-9));18 auto fo = models::fevd_cholesky(fit, 10); for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) CHECK_THAT(fo(i, j), WithinAbs(v["fevd10"][i][j].get<double>(), 1e-9));19 auto g = models::gfevd(fit, 10); for (int i = 0; i < N; ++i) { double rs = 0; for (int j = 0; j < N; ++j) { rs += g(i, j); CHECK_THAT(g(i, j), WithinAbs(v["gfevd10"][i][j].get<double>(), 1e-9)); } CHECK_THAT(rs, WithinAbs(1.0, 1e-12)); }20 auto c = models::connectedness(fit, 10); CHECK_THAT(c.total, WithinAbs(v["dy_total10"].get<double>(), 1e-7)); for (int i = 0; i < N; ++i) { CHECK_THAT(c.from[i], WithinAbs(v["dy_from10"][i].get<double>(), 1e-7)); CHECK_THAT(c.to[i], WithinAbs(v["dy_to10"][i].get<double>(), 1e-7)); }21 auto ord = models::var_select_order(Y, 4); CHECK(ord.bic_best == v["bic_order"].get<int>());22 auto lw = models::ledoit_wolf(Y); CHECK_THAT(lw.shrinkage, WithinAbs(v["lw_shrinkage"].get<double>(), 1e-9)); for (int i = 0; i < N; ++i) CHECK_THAT(lw.cov(i, i), WithinRel(v["lw_cov"][i][i].get<double>(), 1e-9));23}24TEST_CASE("Rolling connectedness runs in parallel and totals are bounded", "[var]") {25 int N; auto Y = load_var_csv(N); auto r = models::rolling_connectedness(Y, 2, 10, 300, 25); CHECK(r.n_windows > 5); CHECK(r.failures == 0); for (double t : r.total) { CHECK(t >= 0); CHECK(t <= 100); }26}27