// E3 — GARCH-family panel estimation on real daily close-to-close returns (2024, from the 1-minute exports):
// per-symbol CPU QML fits (GARCH, GJR, EGARCH; Gaussian and Student-t): convergence rates, iterations, time, SE;
// GPU batched likelihood grid (P = 256 parameter vectors per series) used (i) as a profile-likelihood screen and
// (ii) to choose starting values — effect on convergence and time; 1-day-ahead variance forecasts evaluated against
// next-day RV (QLIKE) on the last 40 sessions (strict chronological split). Usage: bench_e3_garch <e0.json> <out.jsonl> [S]
#include <cstdio>
#include <fstream>
#include <nlohmann/json.hpp>
#include "raid/core/manifest.hpp"
#include "raid/core/hash.hpp"
#include "raid/core/thread_pool.hpp"
#include "raid/data/io.hpp"
#include "raid/data/calendar.hpp"
#include "raid/data/returns.hpp"
#include "raid/estimators/realized.hpp"
#include "raid/estimators/har.hpp"
#include "raid/models/garch.hpp"
#include "raid/backend/backend.hpp"
using namespace raid; using json = nlohmann::json;
namespace { f64 now() { return std::chrono::duration<f64>(std::chrono::steady_clock::now().time_since_epoch()).count(); } }

int main(int argc, char** argv) {
  if (argc < 3) { std::fprintf(stderr, "usage: %s e0_universe.json out.jsonl [S]\n", argv[0]); return 2; }
#ifndef RAID_WITH_ARROW
  return 1;
#else
  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);
  const std::string daily_path = argc > 4 ? argv[4] : "";  // e0_daily.json: 10-year daily bars → long training histories
  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]);
  MarketCalendar cal; std::map<std::string, std::vector<f64>> ret, rvd;
  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()); } }
  // replace the 2024-only close-to-close returns by the 2015–2024 daily history when available (the 2024 intraday RV stays the evaluation target)
  usize replaced = 0; if (!daily_path.empty()) { json ed = json::parse(std::ifstream(daily_path)); man.inputs["e0_daily.json"] = Sha256::of_file(daily_path);
    for (auto& ex : ed["exports"]) { auto pb = data::read_bars_parquet(ex["path"].get<std::string>(), {}, Interval::d1); for (auto& [sym, b] : pb.series) { if (!ret.count(sym) || b.size() < 500) continue; std::vector<f64> r; for (usize i = 1; i < b.size(); ++i) if (b.close[i] > 0 && b.close[i - 1] > 0) r.push_back(std::log(b.close[i] / b.close[i - 1]));
        // keep the last `test` 2024 sessions as evaluation window: the daily history already covers 2024, so align by count: the last 40 returns of the daily series should correspond to the last 40 intraday sessions (both end 2024-12-31)
        ret[sym] = r; ++replaced; } } }
  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);
  out << json({{"kind", "load"}, {"experiment", "E3"}, {"symbols", syms.size()}, {"daily_history_replaced", replaced}}).dump() << "\n";
  const usize S = syms.size(); std::fprintf(stderr, "E3: %zu symbols, ~%zu daily returns each\n", S, ret[syms[0]].size());
  const usize test = 40;  // last 40 sessions held out for forecast evaluation
  // ---- (i) per-symbol fits, all models, CPU parallel over symbols
  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}};
  for (auto& sp : specs) { std::vector<json> rows(S); const f64 t0 = now();
    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()}};
        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;
          // 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 RV
          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]][rvd[syms[i]].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]][rvd[syms[i]].size() - test + k - 1]); auto ln = est::forecast_loss(act, naive); row["qlike_rw_rv"] = ln.qlike; } }
        catch (const RaidError& ex) { row["error"] = ex.what(); row["converged"] = false; } rows[i] = row; } });
    usize conv = 0; for (auto& r : rows) { conv += r.value("converged", false); r["kind"] = "fit"; r["experiment"] = "E3"; out << r.dump() << "\n"; }
    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();
    std::fprintf(stderr, "E3 %s: %zu/%zu converged in %.1fs\n", sp.name, conv, S, now() - t0); }
  // ---- (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/time
  { 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();
      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; } }
    backend::Batch b{cat, off}; auto cpuN = backend::make_cpu_backend(0); std::shared_ptr<backend::Backend> gpu;
#ifdef RAID_ENABLE_METAL
    gpu = backend::make_metal_backend();
#endif
    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";
      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])); }
        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"; } }
    out.flush(); std::fprintf(stderr, "E3 grid screen done\n"); }
  man.save(std::string(argv[2]) + ".manifest.json"); return 0;
#endif
}