1// E2 — HAR-RV on the real daily RV panel (from the 2024 1-minute exports): (a) rolling causal forecasts h ∈ {1,5,22}2// with QLIKE/MSE vs a random-walk benchmark (RV_{t+h} forecast = RV_t) on a strict chronological split, per symbol;3// (b) runtime of per-series QR fits (CPU) vs batched sufficient-statistics fits (CPU-MT, Metal) for S ∈ {10,100,1000}.4// Note: 2024 alone gives ≈ 252 daily RVs per symbol — rolling windows are short (120 days); the point is runtime and5// relative accuracy, not a new forecasting result. Usage: bench_e2_har <e0_universe.json> <out.jsonl> [reps]6#include <cstdio>7#include <fstream>8#include <random>9#include <nlohmann/json.hpp>10#include "raid/core/manifest.hpp"11#include "raid/core/hash.hpp"12#include "raid/core/thread_pool.hpp"13#include "raid/data/io.hpp"14#include "raid/data/calendar.hpp"15#include "raid/data/returns.hpp"16#include "raid/estimators/realized.hpp"17#include "raid/estimators/har.hpp"18using namespace raid; using json = nlohmann::json;19namespace { f64 now() { return std::chrono::duration<f64>(std::chrono::steady_clock::now().time_since_epoch()).count(); } }2021int main(int argc, char** argv) {22 if (argc < 3) { std::fprintf(stderr, "usage: %s e0_universe.json out.jsonl [reps]\n", argv[0]); return 2; }23#ifndef RAID_WITH_ARROW24 return 1;25#else26 json e0 = json::parse(std::ifstream(argv[1])); const int reps = argc > 3 ? std::atoi(argv[3]) : 5; std::ofstream out(argv[2], std::ios::app);27 ExperimentManifest man; man.experiment = "E2.har"; man.created_at_utc = now_utc_iso(); man.raid_version = RAID_VERSION; man.host = host_info(); man.inputs["e0_universe.json"] = Sha256::of_file(argv[1]);28 MarketCalendar cal; std::map<std::string, std::vector<f64>> rv, rq;29 for (auto& ex : e0["exports"]) { auto pb = data::read_bars_parquet(ex["path"].get<std::string>()); for (auto& [sym, b] : pb.series) { auto reg = b.filter_session(cal, SessionPolicy::regular); if (reg.n_sessions() < 240) continue; auto ir = data::intraday_returns(reg); est::RealizedOptions ro; ro.plan.backend = Backend::cpu; auto m = est::realized_measures(ir, ro); bool ok = true; for (f64 v : m.rv) if (!(v > 0)) ok = false; if (ok) { rv[sym] = m.rv; rq[sym] = m.rq; } } }30 const auto universe = e0["samples"]["1000"].get<std::vector<std::string>>(); std::vector<std::string> syms; for (auto& s : universe) if (rv.count(s)) syms.push_back(s);31 std::fprintf(stderr, "E2: %zu symbols with >= 240 sessions of RV\n", syms.size()); out << json({{"kind", "load"}, {"experiment", "E2"}, {"symbols", syms.size()}}).dump() << "\n";32 // (a) forecasting accuracy: window 120, step 1, horizons 1/5/22, HAR-RV vs log-HAR vs HARQ vs random walk; first 100 symbols33 const usize nacc = std::min<usize>(100, syms.size());34 for (int h : {1, 5, 22}) { std::vector<json> rows(nacc);35 ThreadPool::global().parallel_for(nacc, [&](usize a, usize b) { for (usize i = a; i < b; ++i) { auto& x = rv[syms[i]]; auto& q = rq[syms[i]]; json r = {{"symbol", syms[i]}, {"h", h}, {"T", x.size()}};36 for (auto& [name, lg, hq] : std::vector<std::tuple<std::string, bool, bool>>{{"har", false, false}, {"loghar", true, false}, {"harq", false, true}}) { est::HarSpec s; s.horizon = h; s.log = lg; s.harq = hq; try { auto roll = est::har_rolling(x, s, 120, 1, hq ? std::span<const f64>(q) : std::span<const f64>()); auto l = est::forecast_loss(roll.actual, roll.forecast); r[name] = {{"qlike", l.qlike}, {"mse", l.mse}, {"n", l.n}, {"fits", roll.n_fits}, {"seconds", roll.seconds}}; if (name == "har") { std::vector<f64> rw(roll.origin.size()); for (usize j = 0; j < rw.size(); ++j) rw[j] = x[roll.origin[j]]; auto lr = est::forecast_loss(roll.actual, rw); r["random_walk"] = {{"qlike", lr.qlike}, {"mse", lr.mse}, {"n", lr.n}}; } } catch (const RaidError& ex) { r[name] = {{"error", ex.what()}}; } }37 rows[i] = r; } });38 for (auto& r : rows) { r["kind"] = "accuracy"; r["experiment"] = "E2"; out << r.dump() << "\n"; } out.flush(); std::fprintf(stderr, "E2 accuracy h=%d done\n", h); }39 // (b) runtime: per-series fits vs batched sufficient statistics40 for (usize S0 : {10ul, 100ul, 1000ul}) { const usize S = std::min(S0, syms.size()); if (S0 > syms.size() && S0 / 10 >= syms.size()) continue; std::vector<f64> cat; std::vector<i64> off = {0}; for (usize i = 0; i < S; ++i) { cat.insert(cat.end(), rv[syms[i]].begin(), rv[syms[i]].end()); off.push_back(static_cast<i64>(cat.size())); }41 est::HarSpec spec; spec.horizon = 1; spec.cov.type = stats::CovType::classical;42 for (auto mode : {"per_series_cpu1", "per_series_cpuMT", "batched_cpu", "batched_metal"}) for (int rep = 0; rep < reps; ++rep) { const f64 t0 = now(); usize okc = 0; std::string backend_used = "cpu"; f64 dev = 0;43 if (std::string(mode).rfind("per_series", 0) == 0) { std::vector<int> ok(S, 0); auto fitall = [&](usize a, usize b) { for (usize i = a; i < b; ++i) { try { auto e = est::har_fit(rv[syms[i]], spec); ok[i] = e.ok(); } catch (...) {} } }; if (std::string(mode) == "per_series_cpu1") fitall(0, S); else ThreadPool::global().parallel_for(S, fitall); for (int v : ok) okc += v; }44 else { ExecutionPlan p; p.backend = std::string(mode) == "batched_metal" ? Backend::metal : Backend::cpu; auto b = est::har_batch(cat, off, spec, p); for (auto& e : b.fits) okc += e.ok(); backend_used = std::string(to_string(b.stats.backend)); dev = b.stats.device_seconds; }45 out << json({{"kind", "runtime"}, {"experiment", "E2"}, {"S", S}, {"mode", mode}, {"backend_used", backend_used}, {"rep", rep}, {"wall_seconds", now() - t0}, {"device_seconds", dev}, {"converged", okc}, {"T_per_series", cat.size() / S}}).dump() << "\n"; }46 out.flush(); std::fprintf(stderr, "E2 runtime S=%zu done\n", S); }47 man.save(std::string(argv[2]) + ".manifest.json"); return 0;48#endif49}50