1// Monte Carlo evaluation (charter §7): bias, RMSE, empirical coverage of 95 % intervals, convergence/failure rates.2// Sizes are kept small enough for CI (seconds); the full study with more replications runs in E5 (Phase 6).3#include <catch2/catch_test_macros.hpp>4#include <catch2/matchers/catch_matchers_floating_point.hpp>5#include "raid/models/garch.hpp"6#include "raid/estimators/har.hpp"7#include "raid/estimators/ols.hpp"8#include "raid/estimators/realized.hpp"9#include "raid/statistics/bootstrap.hpp"10#include "raid/statistics/reductions.hpp"11#include "raid/core/rng.hpp"12#include "raid/core/thread_pool.hpp"13#include <cmath>14#include <atomic>15using namespace raid; using Catch::Matchers::WithinAbs;1617TEST_CASE("GARCH(1,1) QML: bias, RMSE, coverage of QMLE 95% intervals (R = 120, n = 2000)", "[simulation][garch]") {18 const int R = 120; const usize n = 2000; const double mu = 0.0, omega = 0.02, alpha = 0.08, beta = 0.90;19 std::vector<std::vector<double>> est(R, std::vector<double>(4)); std::vector<std::vector<double>> se(R, std::vector<double>(4)); std::atomic<int> failures{0};20 ThreadPool::global().parallel_for(R, [&](usize a, usize b) { for (usize rep = a; rep < b; ++rep) { Rng rng(777, rep); std::vector<double> r(n); double s2 = omega / (1 - alpha - beta), e = 0;21 for (usize t = 0; t < n; ++t) { s2 = omega + alpha * e * e + beta * s2; e = std::sqrt(s2) * rng.normal(); r[t] = mu + e; }22 models::GarchFitOptions o; auto f = models::garch_fit(r, o); if (!f.ok()) { ++failures; est[rep][0] = NAN; continue; } for (int j = 0; j < 4; ++j) { est[rep][j] = f.coefficients[j]; se[rep][j] = f.inference->std_errors[j]; } } });23 const double truth[4] = {mu, omega, alpha, beta}; const char* names[4] = {"mu", "omega", "alpha", "beta"};24 int ok = 0; for (int rep = 0; rep < R; ++rep) if (std::isfinite(est[rep][0])) ++ok;25 CHECK(failures.load() <= R / 20);26 for (int j = 0; j < 4; ++j) { double bias = 0, mse = 0; int cover = 0, cnt = 0; for (int rep = 0; rep < R; ++rep) { if (!std::isfinite(est[rep][0])) continue; const double d = est[rep][j] - truth[j]; bias += d; mse += d * d; if (std::abs(d) <= 1.96 * se[rep][j]) ++cover; ++cnt; }27 bias /= cnt; const double rmse = std::sqrt(mse / cnt), cov = static_cast<double>(cover) / cnt;28 std::fprintf(stderr, "[mc] garch %-5s bias=%+.4f rmse=%.4f coverage95=%.3f (n=%zu, R=%d, failures=%d)\n", names[j], bias, rmse, cov, n, R, failures.load());29 INFO(names[j]); CHECK(std::abs(bias) < 2.5 * rmse / std::sqrt(static_cast<double>(cnt)) + 0.02); // no bias beyond MC noise + 0.02 (β is known to be slightly downward biased in finite samples)30 CHECK(cov >= 0.85); CHECK(cov <= 0.995); }31}32TEST_CASE("HAR-RV with HAC: coverage of 95% intervals under an AR(1) log-RV DGP with overlapping h = 5 targets (R = 200)", "[simulation][har]") {33 const int R = 200; const usize T = 1200;34 // DGP: log RV_t = c + φ log RV_{t-1} + σ ε_t. The HAR population coefficients are not analytical for an AR(1) in logs;35 // we instead check coverage of the *pseudo-true* value obtained from one very long sample (T = 200 000), which is the36 // relevant target for inference on a misspecified-but-well-defined projection.37 auto simulate = [](usize n, u64 seed) { Rng g(4242, seed); std::vector<double> rv(n); double lv = std::log(1e-4); for (usize t = 0; t < n; ++t) { lv = 0.3 * std::log(1e-4) + 0.7 * lv + 0.4 * g.normal(); rv[t] = std::exp(lv); } return rv; };38 est::HarSpec spec; spec.horizon = 5; auto big = simulate(200000, 999); auto pseudo = est::har_fit(big, spec);39 std::vector<int> cover(4, 0); std::atomic<int> cnt{0}; std::mutex m;40 ThreadPool::global().parallel_for(R, [&](usize a, usize b) { for (usize rep = a; rep < b; ++rep) { auto rv = simulate(T, rep); auto f = est::har_fit(rv, spec); std::lock_guard<std::mutex> lk(m); ++cnt; for (int j = 0; j < 4; ++j) if (std::abs(f.coefficients[j] - pseudo.coefficients[j]) <= 1.96 * f.inference->std_errors[j]) ++cover[j]; } });41 for (int j = 0; j < 4; ++j) { const double c = static_cast<double>(cover[j]) / cnt; std::fprintf(stderr, "[mc] har h=5 coef %d coverage95(HAC)=%.3f (T=%zu, R=%d)\n", j, c, T, R); CHECK(c >= 0.85); }42}43TEST_CASE("Realized variance converges to integrated variance as M grows; semivariances sum to RV", "[simulation][realized]") {44 const double iv = 2.5e-4; for (usize M : {39ul, 390ul, 3900ul}) { double rmse = 0; const int R = 200; for (int rep = 0; rep < R; ++rep) { Rng g(5, static_cast<u64>(M) * 1000 + rep); std::vector<double> r(M); for (auto& v : r) v = std::sqrt(iv / M) * g.normal(); std::vector<i64> off = {0, static_cast<i64>(M)}; est::RealizedOptions o; o.plan.backend = Backend::cpu; auto m = est::realized_measures_batch(backend::Batch{r, off}, o); rmse += (m.rv[0] - iv) * (m.rv[0] - iv); CHECK_THAT(m.rsv_neg[0] + m.rsv_pos[0], WithinAbs(m.rv[0], 1e-18)); }45 rmse = std::sqrt(rmse / R); const double theory = iv * std::sqrt(2.0 / M); std::fprintf(stderr, "[mc] RV rmse M=%zu: %.3e (theory √(2/M)·IV = %.3e)\n", M, rmse, theory); CHECK_THAT(rmse, WithinAbs(theory, 0.3 * theory)); }46}47TEST_CASE("Numerical gradient of the GARCH NLL is stable across step sizes (finite-difference check)", "[simulation][garch]") {48 Rng g(3); const usize n = 1500; std::vector<double> r(n); double s2 = 0.2, e = 0; for (usize t = 0; t < n; ++t) { s2 = 0.02 + 0.08 * e * e + 0.9 * s2; e = std::sqrt(s2) * g.normal(); r[t] = e; }49 std::vector<double> th = {0.0, 0.02, 0.08, 0.9}; auto f = [&](const std::vector<double>& t) { return models::garch_nll(r, t, models::GarchModel::garch, models::GarchDist::gaussian); };50 for (usize j = 0; j < 4; ++j) { double gs[3]; int i = 0; for (double h : {1e-4, 1e-5, 1e-6}) { auto a = th, b = th; a[j] += h; b[j] -= h; gs[i++] = (f(a) - f(b)) / (2 * h); } INFO("param " << j << " grads " << gs[0] << " " << gs[1] << " " << gs[2]); CHECK(std::abs(gs[0] - gs[1]) <= 1e-3 * std::max(1.0, std::abs(gs[1]))); CHECK(std::abs(gs[1] - gs[2]) <= 1e-2 * std::max(1.0, std::abs(gs[1]))); }51}52TEST_CASE("Stationary bootstrap CI for the mean of an AR(1): empirical coverage (R = 100)", "[simulation][bootstrap]") {53 const int R = 100; std::atomic<int> cover{0};54 ThreadPool::global().parallel_for(R, [&](usize a, usize b) { for (usize rep = a; rep < b; ++rep) { Rng g(12, rep); const usize n = 800; std::vector<double> x(n); x[0] = 0; for (usize t = 1; t < n; ++t) x[t] = 0.5 * x[t - 1] + g.normal();55 stats::BootstrapOptions o; o.replicas = 299; o.mean_block = 10; o.seed = 100 + rep; o.threads = 1; auto bs = stats::bootstrap(x, [](std::span<const double> v) { return stats::mean<double>(v); }, o); if (bs.ci_lo <= 0.0 && 0.0 <= bs.ci_hi) ++cover; } });56 const double c = cover.load() / static_cast<double>(R); std::fprintf(stderr, "[mc] stationary bootstrap mean AR(1) coverage95=%.2f\n", c); CHECK(c >= 0.82);57}58