WP40 RAID: Phase 4 — HAR family (fit/rolling/batched), block bootstrap, R references (highfrequency, rugarch), statsmodels HAR refs, Monte Carlo tests, calendar one-offs, Python/CLI HAR+bootstrap
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
22 changed files +2,956 −35
modified
RESEARCH_STATUS.md
+6 −6
@@ -1,6 +1,6 @@ | ||
| 1 | 1 | # RESEARCH_STATUS.md — RAID (UQO WP 40) |
| 2 | 2 | |
| 3 | −Last update: 2026-11-02 (Phase 3 passed). Legend: ✅ verified by tests/artifacts · 🧪 implemented, validation pending · | |
| 3 | +Last update: 2026-11-02 (Phase 4 passed). Legend: ✅ verified by tests/artifacts · 🧪 implemented, validation pending · | |
| 4 | 4 | 🧭 experimental · ⛔ blocked · ⏳ planned · ✖ not implemented (deliberately). |
| 5 | 5 | |
| 6 | 6 | ## Phases |
@@ -10,7 +10,7 @@ Last update: 2026-11-02 (Phase 3 passed). Legend: ✅ verified by tests/artifact | ||
| 10 | 10 | | 1 Foundation | clean build + tests + docs | ✅ passed 2026-11-02 | `ctest` 26/26 (314 assertions), `pytest` 7/7, `docs/PROGRESS.md` Phase 1 | |
| 11 | 11 | | 2 Data | offline synthetic tests + opt-in live tests | ✅ passed 2026-11-02 | offline `[findata]` 4/4, live `[live]` 2/2 (60 assertions), `docs/operations/findata.md` | |
| 12 | 12 | | 3 Numerical primitives | accuracy + measured break-even | ✅ passed 2026-11-02 | `[metal]` 6 cases / 6 411 assertions; `experiments/results/bench_primitives_summary.json` (Metal wins: fused realized ≥ 39 k elements up to 3.4×, GARCH NLL P=64 up to 7.9×; CPU wins: rolling sums, GARCH P=1, small panels) | |
| 13 | −| 4 Econometrics (Tier A) | references + Monte Carlo | ⏳ | — | | |
| 13 | +| 4 Econometrics (Tier A) | references + Monte Carlo | ✅ passed 2026-11-02 | R (highfrequency/rugarch) + statsmodels/arch references, `tests/simulation` Monte Carlo, `docs/methods/registry.md` | | |
| 14 | 14 | | 5 Extensions (Tier B) | per-feature specs + tests | ⏳ | — | |
| 15 | 15 | | 6 Benchmarking E1–E7 | every claim traceable | ⏳ | — | |
| 16 | 16 | | 7 Manuscript | reproducible build, no placeholders | ⏳ | — | |
@@ -22,9 +22,9 @@ Last update: 2026-11-02 (Phase 3 passed). Legend: ✅ verified by tests/artifact | ||
| 22 | 22 | | returns (log/arith), lag/lead, rolling windows, calendar-aware aggregation | A | ✅ | ⏳ (Phase 3) | numpy; exchange_calendars (calendar exact 2000–2026) | ✅ CPU | |
| 23 | 23 | | RV, BPV, RSV±, RQ, jump proxy, TriPV, rCov/rCorr | A | ✅ (numpy truth 1e−12) | ✅ fused kernel 1e−6 vs FP64 | highfrequency 1.0.3 (R comparison Phase 4) | ✅ CPU | |
| 24 | 24 | | OLS/WLS/ridge, HC0–3, HAC, cluster, rolling OLS | A | ✅ (statsmodels 1e−9) | ✅ sufficient statistics 1e−6 (k ≤ 8) | statsmodels 0.15, scipy 1.18 | ✅ CPU | |
| 25 | −| HAR-RV, log-HAR, HAR-X, HARQ | A | ⏳ | ⏳ | highfrequency `HARmodel`, statsmodels | ⏳ | | |
| 25 | +| HAR-RV, log-HAR, HAR-X, HARQ, rolling, batched | A | ✅ (highfrequency 1e−8, statsmodels 1e−9) | ✅ batched sufficient statistics | highfrequency 1.0.3, statsmodels 0.15 | ✅ | | |
| 26 | 26 | | GARCH(1,1), GJR, EGARCH (N, t), QMLE/Hessian/OPG SE, forecasts | A | ✅ (arch 8.0: loglik 5e−4, params 5 %) | ✅ batched NLL 2e−5 (grids) | arch 8.0 ✅, rugarch 1.5.6 (Phase 4) | ✅ CPU | |
| 27 | −| Bootstrap (stationary / circular block), Monte Carlo | A | ⏳ | ⏳ | arch.bootstrap | ⏳ | | |
| 27 | +| Bootstrap (stationary / circular / iid), Monte Carlo harness | A | ✅ (theory + MC coverage) | ⏳ GPU replicas (E5) | theory (Politis–Romano) | ✅ | | |
| 28 | 28 | | VAR, GFEVD, Diebold–Yilmaz | B | ⏳ | ⏳ | statsmodels VAR | ⏳ | |
| 29 | 29 | | Realized GARCH, HEAVY, DCC, PCA/shrinkage, Kalman | B | ✖ until Tier A passes | | | ⏳ | |
| 30 | 30 | | Pre-averaging, realized kernels, MSRV, BNS test | C | ✖ (1-min bars) | | highfrequency | 🧭 only if validated | |
@@ -35,5 +35,5 @@ Last update: 2026-11-02 (Phase 3 passed). Legend: ✅ verified by tests/artifact | ||
| 35 | 35 | not tested unless they come back. |
| 36 | 36 | |
| 37 | 37 | ## Next safe action |
| 38 | −Phase 4: HAR family (fit, rolling, batched via GPU sufficient statistics), block bootstrap, R references | |
| 39 | −(highfrequency / rugarch), Monte Carlo (bias/RMSE/coverage), gradient checks, registry completion. | |
| 38 | +Phase 5 (Tier B, selected for research value): VAR(p) + generalized FEVD + Diebold–Yilmaz connectedness (rolling, batched | |
| 39 | +realized-covariance inputs) for E4; shrinkage covariance; then Phase 6 experiments E1–E7. | |
modified
bindings/python/raid_py.cpp
+29 −0
@@ -12,6 +12,9 @@ | ||
| 12 | 12 | #include "raid/estimators/ols.hpp" |
| 13 | 13 | #include "raid/estimators/realized.hpp" |
| 14 | 14 | #include "raid/models/garch.hpp" |
| 15 | +#include "raid/estimators/har.hpp" | |
| 16 | +#include "raid/statistics/bootstrap.hpp" | |
| 17 | +#include "raid/statistics/reductions.hpp" | |
| 15 | 18 | #include "raid/data/bar_series.hpp" |
| 16 | 19 | #include "raid/data/calendar.hpp" |
| 17 | 20 | #include "raid/data/returns.hpp" |
@@ -94,6 +97,32 @@ NB_MODULE(_core, m) { | ||
| 94 | 97 | |
| 95 | 98 | m.def("calendar_sessions", [](const std::string& a, const std::string& b) { MarketCalendar cal; std::vector<std::string> out; for (auto d : cal.sessions(date_from_iso(a), date_from_iso(b))) out.push_back(d.iso()); return out; }); |
| 96 | 99 | |
| 100 | + // ---- HAR family ---- | |
| 101 | + m.def("har_fit", [](F64 rv, std::vector<usize> lags, int horizon, bool log_model, bool harq, std::optional<F64> rq, const std::string& cov, int hac_lags, const std::string& target) { | |
| 102 | + est::HarSpec s; s.lags = lags; s.horizon = horizon; s.log = log_model; s.harq = harq; s.cov.type = stats::covtype_from_string(cov); s.cov.hac_lags = hac_lags; s.target = target == "single" ? est::HarTarget::single : est::HarTarget::average; | |
| 103 | + std::span<const f64> rvs(rv.data(), rv.shape(0)); std::span<const f64> rqs; if (rq) rqs = std::span<const f64>(rq->data(), rq->shape(0)); | |
| 104 | + json out; { nb::gil_scoped_release rel; auto e = est::har_fit(rvs, s, rqs); out = e.to_json(); out["forecast_next"] = est::har_forecast(e, rvs, s, rqs); } | |
| 105 | + return json_to_py(out); | |
| 106 | + }, nb::arg("rv"), nb::arg("lags") = std::vector<usize>{1, 5, 22}, nb::arg("horizon") = 1, nb::arg("log") = false, nb::arg("harq") = false, nb::arg("rq") = nb::none(), nb::arg("cov") = "hac", nb::arg("hac_lags") = -1, nb::arg("target") = "average"); | |
| 107 | + m.def("har_rolling", [](F64 rv, std::vector<usize> lags, int horizon, bool log_model, usize window, usize step, bool harq, std::optional<F64> rq) { | |
| 108 | + est::HarSpec s; s.lags = lags; s.horizon = horizon; s.log = log_model; s.harq = harq; std::span<const f64> rvs(rv.data(), rv.shape(0)); std::span<const f64> rqs; if (rq) rqs = std::span<const f64>(rq->data(), rq->shape(0)); | |
| 109 | + est::HarRolling r; { nb::gil_scoped_release rel; r = est::har_rolling(rvs, s, window, step, rqs); } | |
| 110 | + nb::dict d; std::vector<int64_t> o(r.origin.begin(), r.origin.end()); d["origin"] = to_np_i64(std::move(o)); d["forecast"] = to_np(r.forecast); d["actual"] = to_np(r.actual); d["n_fits"] = r.n_fits; d["seconds"] = r.seconds; | |
| 111 | + auto l = est::forecast_loss(r.actual, r.forecast); d["mse"] = l.mse; d["qlike"] = l.qlike; d["n_loss"] = l.n; return d; | |
| 112 | + }, nb::arg("rv"), nb::arg("lags") = std::vector<usize>{1, 5, 22}, nb::arg("horizon") = 1, nb::arg("log") = false, nb::arg("window") = 1000, nb::arg("step") = 1, nb::arg("harq") = false, nb::arg("rq") = nb::none()); | |
| 113 | + m.def("har_batch", [](F64 rv, I64 off, std::vector<usize> lags, int horizon, bool log_model, const std::string& backend_s) { | |
| 114 | + est::HarSpec s; s.lags = lags; s.horizon = horizon; s.log = log_model; ExecutionPlan p; p.backend = backend_from_string(backend_s); | |
| 115 | + est::HarBatchResult r; { nb::gil_scoped_release rel; r = est::har_batch(std::span<const f64>(rv.data(), rv.shape(0)), std::span<const i64>(off.data(), off.shape(0)), s, p); } | |
| 116 | + nb::list fits; for (auto& e : r.fits) fits.append(json_to_py(e.to_json())); nb::dict d; d["fits"] = fits; d["backend"] = std::string(to_string(r.stats.backend)); d["precision"] = std::string(to_string(r.stats.precision)); d["wall_seconds"] = r.stats.wall_seconds; d["device_seconds"] = r.stats.device_seconds; return d; | |
| 117 | + }, nb::arg("rv"), nb::arg("offsets"), nb::arg("lags") = std::vector<usize>{1, 5, 22}, nb::arg("horizon") = 1, nb::arg("log") = false, nb::arg("backend") = "auto"); | |
| 118 | + // ---- bootstrap ---- | |
| 119 | + m.def("bootstrap", [](F64 x, const std::string& stat, const std::string& scheme, double mean_block, usize replicas, uint64_t seed, double level) { | |
| 120 | + stats::BootstrapOptions o; o.scheme = scheme == "circular" ? stats::BootstrapScheme::circular_block : (scheme == "iid" ? stats::BootstrapScheme::iid : stats::BootstrapScheme::stationary); o.mean_block = mean_block; o.replicas = replicas; o.seed = seed; o.ci_level = level; | |
| 121 | + std::function<f64(std::span<const f64>)> f; if (stat == "mean") f = [](std::span<const f64> v) { return stats::mean<f64>(v); }; else if (stat == "rv" || stat == "sum_squares") f = [](std::span<const f64> v) { return stats::sum_squares<f64>(v); }; else if (stat == "variance") f = [](std::span<const f64> v) { return stats::variance<f64>(v); }; else throw RaidError(ErrorCode::invalid_argument, "stat must be mean|rv|variance"); | |
| 122 | + stats::BootstrapResult r; { nb::gil_scoped_release rel; r = stats::bootstrap(std::span<const f64>(x.data(), x.shape(0)), f, o); } | |
| 123 | + nb::dict d; d["point"] = r.point; d["se"] = r.se; d["ci_lo"] = r.ci_lo; d["ci_hi"] = r.ci_hi; d["bias"] = r.bias; d["failures"] = r.failures; d["block_length"] = r.block_length; d["statistics"] = to_np(r.statistics); d["seconds"] = r.seconds; return d; | |
| 124 | + }, nb::arg("x"), nb::arg("stat") = "mean", nb::arg("scheme") = "stationary", nb::arg("mean_block") = 0.0, nb::arg("replicas") = 999, nb::arg("seed") = 20261102, nb::arg("level") = 0.95); | |
| 125 | + | |
| 97 | 126 | // ---- FinData client (keys from FINDATA_API_KEY / FINDATA_API_KEY_FILE; never passed through Python) ---- |
| 98 | 127 | m.def("findata_bars", [](const std::string& symbol, const std::string& interval, const std::string& start, const std::string& end, const std::string& adjustment, const std::string& cache_dir, bool offline) { |
| 99 | 128 | findata::Config cfg = findata::Config::from_env(); if (!cache_dir.empty()) cfg.cache_dir = cache_dir; cfg.offline = offline; findata::Client cl(cfg); |
modified
cli/raid.cpp
+10 −0
@@ -18,6 +18,7 @@ | ||
| 18 | 18 | #include "raid/estimators/realized.hpp" |
| 19 | 19 | #include "raid/estimators/ols.hpp" |
| 20 | 20 | #include "raid/models/garch.hpp" |
| 21 | +#include "raid/estimators/har.hpp" | |
| 21 | 22 | #ifdef RAID_WITH_FINDATA |
| 22 | 23 | #include "raid/data/findata.hpp" |
| 23 | 24 | #endif |
@@ -34,6 +35,7 @@ void usage() { | ||
| 34 | 35 | " raid estimate realized <bars> [--backend auto|cpu|metal] [--bpv-no-fs] per-session realized measures\n" |
| 35 | 36 | " raid estimate garch <bars> [--model garch|gjr|egarch] [--dist gaussian|student_t] [--horizon H] [--returns daily|intraday]\n" |
| 36 | 37 | " raid estimate ols <csv y,x1,x2...> [--cov classical|hc0..hc3|hac|cluster] [--lags L]\n" |
| 38 | + " raid estimate har <bars> [--horizon h] [--log] [--harq] [--window W --step S] daily RV from 1m bars → HAR fit (+ rolling forecasts)\n" | |
| 37 | 39 | " raid data ... FinData access (Phase 2)\n" |
| 38 | 40 | " raid benchmark ... benchmark harness (Phase 3+)\n", RAID_VERSION); |
| 39 | 41 | } |
@@ -89,6 +91,14 @@ int cmd_estimate(const std::vector<std::string>& a) { | ||
| 89 | 91 | try { auto e = models::garch_fit(r, go); out[sym] = e.to_json(); } catch (const RaidError& ex) { out[sym] = {{"error", ex.what()}}; } } |
| 90 | 92 | std::cout << out.dump(1) << "\n"; return 0; |
| 91 | 93 | } |
| 94 | + if (what == "har") { | |
| 95 | + auto bars = load_bars(path); MarketCalendar cal; json out = json::object(); est::HarSpec s; s.horizon = std::stoi(arg(a, "--horizon", "1")); s.log = flag(a, "--log"); s.harq = flag(a, "--harq"); | |
| 96 | + 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 ro; ro.plan = plan; auto m = est::realized_measures(ir, ro); | |
| 97 | + try { auto e = est::har_fit(m.rv, s, m.rq); json j = e.to_json(); j["forecast_next"] = est::har_forecast(e, m.rv, s, m.rq); j["n_sessions"] = m.rv.size(); | |
| 98 | + if (!arg(a, "--window").empty()) { auto r = est::har_rolling(m.rv, s, static_cast<usize>(std::stoul(arg(a, "--window"))), static_cast<usize>(std::stoul(arg(a, "--step", "1"))), m.rq); auto l = est::forecast_loss(r.actual, r.forecast); j["rolling"] = {{"n_fits", r.n_fits}, {"mse", l.mse}, {"qlike", l.qlike}, {"seconds", r.seconds}}; } | |
| 99 | + out[sym] = j; } catch (const RaidError& ex) { out[sym] = {{"error", ex.what()}}; } } | |
| 100 | + std::cout << out.dump(1) << "\n"; return 0; | |
| 101 | + } | |
| 92 | 102 | if (what == "ols") { |
| 93 | 103 | std::ifstream f(path); RAID_REQUIRE(f.good(), io, "cannot open " + path); std::string line; std::getline(f, line); |
| 94 | 104 | std::vector<std::string> names; { usize s = 0; for (usize i = 0; i <= line.size(); ++i) if (i == line.size() || line[i] == ',') { names.push_back(line.substr(s, i - s)); s = i + 1; } } |
added
datasets/fixtures/har_rv.csv
+1501 −0
Large diff (1,501 lines) is not loaded by default to keep the page fast.Load diff
added
datasets/fixtures/reference_r.json
+275 −0
Large diff (275 lines) is not loaded by default to keep the page fast.Load diff
modified
datasets/fixtures/reference_values.json
+676 −9
Large diff (704 lines) is not loaded by default to keep the page fast.Load diff
modified
docs/PROGRESS.md
+29 −0
Large diff (32 lines) is not loaded by default to keep the page fast.Load diff
modified
docs/methods/registry.md
+12 −8
Large diff (35 lines) is not loaded by default to keep the page fast.Load diff
modified
include/raid/estimators/har.hpp
+3 −2
Large diff (22 lines) is not loaded by default to keep the page fast.Load diff
added
include/raid/statistics/bootstrap.hpp
+55 −0
Large diff (55 lines) is not loaded by default to keep the page fast.Load diff
modified
python/raid/__init__.py
+2 −1
Large diff (9 lines) is not loaded by default to keep the page fast.Load diff
modified
python/raid/api.py
+20 −0
Large diff (23 lines) is not loaded by default to keep the page fast.Load diff
modified
python/tests/test_api.py
+23 −0
Large diff (26 lines) is not loaded by default to keep the page fast.Load diff
modified
scripts/make_fixtures.py
+27 −4
Large diff (41 lines) is not loaded by default to keep the page fast.Load diff
added
scripts/reference_r.R
+50 −0
Large diff (50 lines) is not loaded by default to keep the page fast.Load diff
modified
src/estimators/har.cpp
+7 −3
Large diff (37 lines) is not loaded by default to keep the page fast.Load diff
added
src/statistics/bootstrap.cpp
+48 −0
Large diff (48 lines) is not loaded by default to keep the page fast.Load diff
added
tests/reference/test_har_reference.cpp
+37 −0
Large diff (37 lines) is not loaded by default to keep the page fast.Load diff
added
tests/reference/test_r_reference.cpp
+55 −0
Large diff (55 lines) is not loaded by default to keep the page fast.Load diff
modified
tests/reference/test_reference_values.cpp
+2 −2
Large diff (16 lines) is not loaded by default to keep the page fast.Load diff
added
tests/simulation/test_monte_carlo.cpp
+57 −0
Large diff (57 lines) is not loaded by default to keep the page fast.Load diff
added
tests/unit/test_bootstrap.cpp
+32 −0
Large diff (33 lines) is not loaded by default to keep the page fast.Load diff