1#include "raid/estimators/har.hpp"2#include "raid/core/error.hpp"3#include "raid/statistics/covariance.hpp"4#include <chrono>5#include <cmath>6#include <algorithm>78namespace raid::est {910HarDesign har_design(std::span<const f64> rv, const HarSpec& spec, std::span<const f64> rq, const Eigen::MatrixXd& exog, const std::vector<std::string>& exog_names) {11 const usize T = rv.size(); RAID_REQUIRE(!spec.lags.empty() && spec.horizon >= 1, invalid_argument, "har: lags and horizon");12 const usize maxlag = *std::max_element(spec.lags.begin(), spec.lags.end()); const usize h = static_cast<usize>(spec.horizon);13 RAID_REQUIRE(T > maxlag + h + 2, insufficient_samples, "har: need T > maxlag + h + 2");14 for (f64 v : rv) RAID_REQUIRE(std::isfinite(v) && (!spec.log || v > 0), data_quality, "har: RV must be finite (and > 0 for log-HAR)");15 if (spec.harq) RAID_REQUIRE(rq.size() == T, invalid_argument, "harq: RQ series required with the same length");16 const usize nex = exog.size() ? static_cast<usize>(exog.cols()) : 0; if (nex) RAID_REQUIRE(static_cast<usize>(exog.rows()) == T, invalid_argument, "har-x: exog rows must equal T");17 std::vector<f64> z(T); for (usize t = 0; t < T; ++t) z[t] = spec.log ? std::log(rv[t]) : rv[t];18 f64 sqrt_rq_mean = 0.0;19 if (spec.harq && spec.harq_demean_sqrt_rq) { f64 s = 0; for (usize t = 0; t < T; ++t) s += std::max(rq[t], 0.0); sqrt_rq_mean = std::sqrt(s / static_cast<f64>(T)); } // √(mean RQ) over the sample, as highfrequency::HARmodel20 // prefix sums for the window means (FP64, Neumaier)21 std::vector<f64> cs(T + 1, 0.0); { f64 s = 0, c = 0; for (usize t = 0; t < T; ++t) { const f64 v = z[t], tt = s + v; if (std::abs(s) >= std::abs(v)) c += (s - tt) + v; else c += (v - tt) + s; s = tt; cs[t + 1] = s + c; } }22 auto wmean = [&](usize t, usize k) { f64 s = 0; for (usize j = 0; j < k; ++j) s += z[t - j]; return s / static_cast<f64>(k); }; // direct sum (exact enough, k ≤ 22)23 const usize first = maxlag - 1, last = T - 1 - h; // origins t with full lags and a full target24 HarDesign d; d.first_t = first; d.n = last - first + 1; const usize k = 1 + spec.lags.size() + (spec.harq ? 1 : 0) + nex;25 d.X.resize(static_cast<Eigen::Index>(d.n), static_cast<Eigen::Index>(k)); d.y.resize(static_cast<Eigen::Index>(d.n));26 d.names.push_back("const"); for (usize l : spec.lags) d.names.push_back((spec.log ? "logrv_" : "rv_") + std::to_string(l)); if (spec.harq) d.names.push_back(spec.harq_demean_sqrt_rq ? "sqrt_rq_dm_x_rv" : "sqrt_rq_x_rv");27 d.sqrt_rq_mean = sqrt_rq_mean; for (usize j = 0; j < nex; ++j) d.names.push_back(j < exog_names.size() ? exog_names[j] : "x" + std::to_string(j));28 for (usize t = first; t <= last; ++t) {29 const Eigen::Index i = static_cast<Eigen::Index>(t - first); Eigen::Index c = 0; d.X(i, c++) = 1.0;30 for (usize l : spec.lags) d.X(i, c++) = wmean(t, l);31 if (spec.harq) d.X(i, c++) = (std::sqrt(std::max(rq[t], 0.0)) - sqrt_rq_mean) * rv[t];32 for (usize j = 0; j < nex; ++j) d.X(i, c++) = exog(static_cast<Eigen::Index>(t), static_cast<Eigen::Index>(j));33 if (spec.target == HarTarget::single) d.y[i] = z[t + h]; else { f64 s = 0; for (usize j = 1; j <= h; ++j) s += spec.log ? rv[t + j] : z[t + j]; s /= static_cast<f64>(h); d.y[i] = spec.log ? std::log(s) : s; }34 }35 return d;36}3738Estimate har_fit(std::span<const f64> rv, const HarSpec& spec, std::span<const f64> rq, const Eigen::MatrixXd& exog, const std::vector<std::string>& exog_names) {39 auto d = har_design(rv, spec, rq, exog, exog_names);40 OlsOptions o; o.cov = spec.cov; o.names = d.names;41 if (o.cov.type == stats::CovType::hac_bartlett && o.cov.hac_lags < 0) o.cov.hac_lags = std::max(spec.horizon - 1, stats::default_hac_lags(d.n));42 auto r = ols(d.X, d.y, o); Estimate e = r.estimate; e.model = spec.model_key();43 e.extra["horizon"] = spec.horizon; e.extra["target"] = spec.target == HarTarget::average ? "average" : "single"; e.extra["lags"] = spec.lags; e.extra["first_origin"] = d.first_t; e.extra["log"] = spec.log; e.extra["harq"] = spec.harq;44 if (spec.log) { f64 s = 0; for (Eigen::Index i = 0; i < r.residuals.size(); ++i) s += std::exp(r.residuals[i]); e.extra["smearing"] = s / static_cast<f64>(r.residuals.size()); }45 if (spec.harq) e.extra["sqrt_rq_mean"] = d.sqrt_rq_mean;46 return e;47}4849f64 har_forecast(const Estimate& fit, std::span<const f64> rv, const HarSpec& spec, std::span<const f64> rq, const Eigen::VectorXd& exog_last) {50 const usize T = rv.size(); const usize maxlag = *std::max_element(spec.lags.begin(), spec.lags.end()); RAID_REQUIRE(T >= maxlag, insufficient_samples, "har_forecast: not enough history");51 const usize t = T - 1; f64 yhat = fit.coefficients[0]; usize c = 1;52 for (usize l : spec.lags) { f64 s = 0; for (usize j = 0; j < l; ++j) s += spec.log ? std::log(rv[t - j]) : rv[t - j]; yhat += fit.coefficients[c++] * s / static_cast<f64>(l); }53 if (spec.harq) { RAID_REQUIRE(rq.size() == T, invalid_argument, "harq forecast needs RQ"); const f64 dm = fit.extra.value("sqrt_rq_mean", 0.0); yhat += fit.coefficients[c++] * (std::sqrt(std::max(rq[t], 0.0)) - dm) * rv[t]; }54 for (Eigen::Index j = 0; j < exog_last.size(); ++j) yhat += fit.coefficients[c++] * exog_last[j];55 return spec.log ? std::exp(yhat) : yhat;56}5758HarRolling har_rolling(std::span<const f64> rv, const HarSpec& spec, usize window, usize step, std::span<const f64> rq) {59 const auto t0 = std::chrono::steady_clock::now(); HarRolling out; const usize T = rv.size(); const usize h = static_cast<usize>(spec.horizon);60 const usize maxlag = *std::max_element(spec.lags.begin(), spec.lags.end()); const usize minT = maxlag + h + 30;61 RAID_REQUIRE(window == 0 || window >= minT, insufficient_samples, "har_rolling: window too small"); RAID_REQUIRE(T > minT + h, insufficient_samples, "har_rolling: series too short");62 const usize k = 1 + spec.lags.size() + (spec.harq ? 1 : 0); std::vector<Eigen::VectorXd> coefs;63 for (usize end = (window ? window : minT); end + h <= T; end += step) { // model uses rv[0..end) (or the last `window`), forecasts target at origin end−164 const usize b = window ? end - window : 0; std::span<const f64> hist = rv.subspan(b, end - b); std::span<const f64> rqh = spec.harq ? rq.subspan(b, end - b) : std::span<const f64>();65 Estimate e; try { HarSpec s2 = spec; s2.cov.type = stats::CovType::classical; e = har_fit(hist, s2, rqh); } catch (const RaidError&) { continue; }66 const f64 fc = har_forecast(e, hist, spec, rqh); const usize t = end - 1;67 f64 actual; if (spec.target == HarTarget::single) actual = rv[t + h]; else { f64 s = 0; for (usize j = 1; j <= h; ++j) s += rv[t + j]; actual = s / static_cast<f64>(h); }68 out.origin.push_back(t); out.forecast.push_back(fc); out.actual.push_back(actual); coefs.push_back(Eigen::Map<const Eigen::VectorXd>(e.coefficients.data(), static_cast<Eigen::Index>(k))); ++out.n_fits;69 }70 out.coefficients.resize(static_cast<Eigen::Index>(k), static_cast<Eigen::Index>(coefs.size())); for (usize j = 0; j < coefs.size(); ++j) out.coefficients.col(static_cast<Eigen::Index>(j)) = coefs[j];71 out.seconds = std::chrono::duration<f64>(std::chrono::steady_clock::now() - t0).count(); return out;72}7374HarBatchResult har_batch(std::span<const f64> rv_concat, std::span<const i64> offsets, const HarSpec& spec, const ExecutionPlan& plan) {75 HarBatchResult out; const usize S = offsets.size() - 1; std::string why; auto be = backend::resolve(plan, &why);76 // Build all designs (CPU, FP64), concatenate row-major for the backend, accumulate X'X / X'y / y'y per series, solve in FP64.77 std::vector<HarDesign> designs(S); usize k = 0, rows = 0; std::vector<i64> roff(S + 1, 0); std::vector<usize> ok(S, 1);78 for (usize s = 0; s < S; ++s) { try { designs[s] = har_design(rv_concat.subspan(static_cast<usize>(offsets[s]), static_cast<usize>(offsets[s + 1] - offsets[s])), spec); k = static_cast<usize>(designs[s].X.cols()); rows += designs[s].n; } catch (const RaidError&) { ok[s] = 0; } roff[s + 1] = static_cast<i64>(rows); }79 RAID_REQUIRE(k > 0, insufficient_samples, "har_batch: no series with enough data");80 std::vector<f64> X(rows * k), y(rows); for (usize s = 0; s < S; ++s) { if (!ok[s]) continue; const usize b = static_cast<usize>(roff[s]); for (usize i = 0; i < designs[s].n; ++i) { y[b + i] = designs[s].y[static_cast<Eigen::Index>(i)]; for (usize j = 0; j < k; ++j) X[(b + i) * k + j] = designs[s].X(static_cast<Eigen::Index>(i), static_cast<Eigen::Index>(j)); } }81 std::vector<f64> xtx, xty, yty; be->regression_stats(X, y, roff, k, xtx, xty, yty, &out.stats);82 out.fits.resize(S);83 for (usize s = 0; s < S; ++s) {84 if (!ok[s]) { out.fits[s].model = spec.model_key(); out.fits[s].convergence = Convergence{false, 0, 0, 0, 1, "insufficient data"}; continue; }85 SufficientStats ss; ss.xtx = Eigen::Map<const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>(xtx.data() + s * k * k, static_cast<Eigen::Index>(k), static_cast<Eigen::Index>(k)); ss.xty = Eigen::Map<const Eigen::VectorXd>(xty.data() + s * k, static_cast<Eigen::Index>(k)); ss.yty = yty[s]; ss.n = designs[s].n;86 OlsOptions o; o.names = designs[s].names; auto r = ols_from_stats(ss, o); out.fits[s] = r.estimate; out.fits[s].model = spec.model_key(); out.fits[s].backend = be->kind(); out.fits[s].precision = out.stats.precision;87 out.fits[s].extra["horizon"] = spec.horizon; out.fits[s].extra["note"] = "classical SE from sufficient statistics; use har_fit for HAC";88 }89 return out;90}9192ForecastLoss forecast_loss(std::span<const f64> a, std::span<const f64> f) {93 RAID_REQUIRE(a.size() == f.size(), invalid_argument, "forecast_loss: size mismatch"); ForecastLoss l{0, 0, 0};94 for (usize i = 0; i < a.size(); ++i) { if (!(a[i] > 0) || !(f[i] > 0) || !std::isfinite(a[i]) || !std::isfinite(f[i])) continue; const f64 d = a[i] - f[i]; l.mse += d * d; const f64 q = a[i] / f[i]; l.qlike += q - std::log(q) - 1.0; ++l.n; }95 if (l.n) { l.mse /= static_cast<f64>(l.n); l.qlike /= static_cast<f64>(l.n); } return l;96}9798} // namespace raid::est99