1// E3 — GARCH-family panel estimation on real daily close-to-close returns (2024, from the 1-minute exports):2// per-symbol CPU QML fits (GARCH, GJR, EGARCH; Gaussian and Student-t): convergence rates, iterations, time, SE;3// GPU batched likelihood grid (P = 256 parameter vectors per series) used (i) as a profile-likelihood screen and4// (ii) to choose starting values — effect on convergence and time; 1-day-ahead variance forecasts evaluated against5// next-day RV (QLIKE) on the last 40 sessions (strict chronological split). Usage: bench_e3_garch <e0.json> <out.jsonl> [S]6#include <cstdio>7#include <fstream>8#include <nlohmann/json.hpp>9#include "raid/core/manifest.hpp"10#include "raid/core/hash.hpp"11#include "raid/core/thread_pool.hpp"12#include "raid/data/io.hpp"13#include "raid/data/calendar.hpp"14#include "raid/data/returns.hpp"15#include "raid/estimators/realized.hpp"16#include "raid/estimators/har.hpp"17#include "raid/models/garch.hpp"18#include "raid/backend/backend.hpp"19using namespace raid; using json = nlohmann::json;20namespace { f64 now() { return std::chrono::duration<f64>(std::chrono::steady_clock::now().time_since_epoch()).count(); } }2122int main(int argc, char** argv) {23 if (argc < 3) { std::fprintf(stderr, "usage: %s e0_universe.json out.jsonl [S]\n", argv[0]); return 2; }24#ifndef RAID_WITH_ARROW25 return 1;26#else27 json e0 = json::parse(std::ifstream(argv[1])); const usize Smax = argc > 3 ? static_cast<usize>(std::atoi(argv[3])) : 1000; std::ofstream out(argv[2], std::ios::app);28 ExperimentManifest man; man.experiment = "E3.garch"; 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]);29 MarketCalendar cal; std::map<std::string, std::vector<f64>> ret, rvd;30 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); std::vector<f64> r; for (usize k = 1; k < ir.close_to_close.size(); ++k) r.push_back(ir.close_to_close[k]); est::RealizedOptions ro; ro.plan.backend = Backend::cpu; auto m = est::realized_measures(ir, ro); ret[sym] = r; rvd[sym] = std::vector<f64>(m.rv.begin() + 1, m.rv.end()); } }31 const auto universe = e0["samples"]["1000"].get<std::vector<std::string>>(); std::vector<std::string> syms; for (auto& s : universe) if (ret.count(s) && syms.size() < Smax) syms.push_back(s);32 const usize S = syms.size(); std::fprintf(stderr, "E3: %zu symbols, ~%zu daily returns each\n", S, ret[syms[0]].size());33 const usize test = 40; // last 40 sessions held out for forecast evaluation34 // ---- (i) per-symbol fits, all models, CPU parallel over symbols35 struct Spec { const char* name; models::GarchModel m; models::GarchDist d; }; const std::vector<Spec> specs = {{"garch.n", models::GarchModel::garch, models::GarchDist::gaussian}, {"gjr.n", models::GarchModel::gjr, models::GarchDist::gaussian}, {"egarch.n", models::GarchModel::egarch, models::GarchDist::gaussian}, {"garch.t", models::GarchModel::garch, models::GarchDist::student_t}};36 for (auto& sp : specs) { std::vector<json> rows(S); const f64 t0 = now();37 ThreadPool::global().parallel_for(S, [&](usize a, usize b) { for (usize i = a; i < b; ++i) { auto& r = ret[syms[i]]; std::vector<f64> train(r.begin(), r.end() - static_cast<long>(test)); models::GarchFitOptions o; o.model = sp.m; o.dist = sp.d; o.forecast_horizon = 1; json row = {{"symbol", syms[i]}, {"model", sp.name}, {"n_train", train.size()}};38 try { const f64 f0 = now(); auto e = models::garch_fit(train, o); row["seconds"] = now() - f0; row["converged"] = e.ok(); row["iterations"] = e.convergence.iterations; row["fevals"] = e.convergence.function_evaluations; row["loglik"] = -e.loss; row["coef"] = e.coefficients; row["se"] = e.inference ? e.inference->std_errors : std::vector<f64>{}; row["persistence"] = e.extra.value("persistence", 0.0); row["warnings"] = e.warnings;39 // rolling 1-step forecasts over the test window WITHOUT refitting (parameters fixed at the training estimate; variance path updated with realized returns): QLIKE vs next-day RV40 if (e.ok()) { std::vector<f64> fc, act; for (usize k = 0; k < test; ++k) { std::vector<f64> hist(r.begin(), r.end() - static_cast<long>(test - k)); auto f = models::garch_forecast(e, hist, 1); fc.push_back(f[0]); act.push_back(rvd[syms[i]][r.size() - test + k]); } auto l = est::forecast_loss(act, fc); row["qlike_vs_rv"] = l.qlike; row["mse_vs_rv"] = l.mse; row["n_forecasts"] = l.n; std::vector<f64> naive; for (usize k = 0; k < test; ++k) naive.push_back(rvd[syms[i]][r.size() - test + k - 1]); auto ln = est::forecast_loss(act, naive); row["qlike_rw_rv"] = ln.qlike; } }41 catch (const RaidError& ex) { row["error"] = ex.what(); row["converged"] = false; } rows[i] = row; } });42 usize conv = 0; for (auto& r : rows) { conv += r.value("converged", false); r["kind"] = "fit"; r["experiment"] = "E3"; out << r.dump() << "\n"; }43 out << json({{"kind", "fit_summary"}, {"experiment", "E3"}, {"model", sp.name}, {"S", S}, {"converged", conv}, {"wall_seconds_all", now() - t0}, {"threads", ThreadPool::global().size()}}).dump() << "\n"; out.flush();44 std::fprintf(stderr, "E3 %s: %zu/%zu converged in %.1fs\n", sp.name, conv, S, now() - t0); }45 // ---- (ii) GPU grid screen: P = 256 (α, β, ω) combinations per series (μ = sample mean), batched NLL on Metal vs CPU-MT; best grid point as starting value → refit convergence/time46 { const usize P = 256; std::vector<f64> cat; std::vector<i64> off = {0}; std::vector<f64> params(S * P * 4); for (usize i = 0; i < S; ++i) { auto& r = ret[syms[i]]; std::vector<f64> train(r.begin(), r.end() - static_cast<long>(test)); cat.insert(cat.end(), train.begin(), train.end()); off.push_back(static_cast<i64>(cat.size())); f64 mu = 0, v = 0; for (f64 x : train) mu += x; mu /= train.size(); for (f64 x : train) v += (x - mu) * (x - mu); v /= train.size();47 usize q = 0; for (int ia = 0; ia < 16; ++ia) for (int ib = 0; ib < 16; ++ib) { const f64 alpha = 0.01 + 0.015 * ia, beta = 0.50 + 0.031 * ib; f64* th = params.data() + (i * P + q) * 4; th[0] = mu; th[2] = alpha; th[3] = beta; th[1] = alpha + beta < 0.999 ? v * (1 - alpha - beta) : v * 0.001; ++q; } }48 backend::Batch b{cat, off}; auto cpuN = backend::make_cpu_backend(0); std::shared_ptr<backend::Backend> gpu;49#ifdef RAID_ENABLE_METAL50 gpu = backend::make_metal_backend();51#endif52 std::vector<f64> nll_cpu, nll_gpu; for (int rep = 0; rep < 3; ++rep) { backend::KernelStats st; const f64 t0 = now(); nll_cpu = cpuN->garch_nll_batch(b, params, P, 4, 0, 0.0, &st); out << json({{"kind", "grid"}, {"experiment", "E3"}, {"backend", "cpuMT"}, {"S", S}, {"P", P}, {"rep", rep}, {"wall_seconds", now() - t0}}).dump() << "\n";53 if (gpu) { const f64 t1 = now(); nll_gpu = gpu->garch_nll_batch(b, params, P, 4, 0, 0.0, &st); f64 maxrel = 0; usize agree = 0; for (usize i = 0; i < S; ++i) { usize ac = 0, ag = 0; for (usize q = 1; q < P; ++q) { if (nll_cpu[i * P + q] < nll_cpu[i * P + ac]) ac = q; if (nll_gpu[i * P + q] < nll_gpu[i * P + ag]) ag = q; } agree += (ac == ag); for (usize q = 0; q < P; ++q) if (std::isfinite(nll_cpu[i * P + q]) && std::isfinite(nll_gpu[i * P + q])) maxrel = std::max(maxrel, std::abs(nll_cpu[i * P + q] - nll_gpu[i * P + q]) / std::abs(nll_cpu[i * P + q])); }54 out << json({{"kind", "grid"}, {"experiment", "E3"}, {"backend", "metal"}, {"S", S}, {"P", P}, {"rep", rep}, {"wall_seconds", now() - t1}, {"device_seconds", st.device_seconds}, {"max_rel_err_nll", maxrel}, {"argmin_agreement", agree}, {"S_total", S}}).dump() << "\n"; } }55 out.flush(); std::fprintf(stderr, "E3 grid screen done\n"); }56 man.save(std::string(argv[2]) + ".manifest.json"); return 0;57#endif58}59