1#include "raid/models/garch.hpp"2#include "raid/core/error.hpp"3#include "raid/statistics/reductions.hpp"4#include <cmath>5#include <limits>6#include <numbers>78namespace raid::models {910constexpr f64 kNuMin = 2.05, kNuMax = 500.0; // same bounds as the arch package (StudentsT distribution bounds)1112std::string to_string(GarchModel m) { switch (m) { case GarchModel::garch: return "garch"; case GarchModel::gjr: return "gjr"; case GarchModel::egarch: return "egarch"; } return "?"; }13std::string to_string(GarchDist d) { return d == GarchDist::gaussian ? "gaussian" : "student_t"; }1415usize garch_n_params(GarchModel m, GarchDist d) { return 4 + (m != GarchModel::garch ? 1 : 0) + (d == GarchDist::student_t ? 1 : 0); }16std::vector<std::string> garch_param_names(GarchModel m, GarchDist d) {17 std::vector<std::string> n = {"mu", "omega", "alpha"}; if (m != GarchModel::garch) n.push_back("gamma"); n.push_back("beta"); if (d == GarchDist::student_t) n.push_back("nu"); return n;18}1920bool garch_params_valid(std::span<const f64> th, GarchModel m, GarchDist d) {21 if (th.size() != garch_n_params(m, d)) return false;22 for (f64 v : th) if (!std::isfinite(v)) return false;23 const f64 omega = th[1], alpha = th[2];24 const f64 gamma = m != GarchModel::garch ? th[3] : 0.0;25 const f64 beta = th[m != GarchModel::garch ? 4 : 3];26 if (d == GarchDist::student_t && th.back() <= 2.0) return false;27 switch (m) {28 case GarchModel::garch: return omega > 0 && alpha >= 0 && beta >= 0 && alpha + beta < 1.0;29 case GarchModel::gjr: return omega > 0 && alpha >= 0 && beta >= 0 && alpha + gamma >= 0 && alpha + beta + 0.5 * gamma < 1.0;30 case GarchModel::egarch: return std::abs(beta) < 1.0;31 }32 return false;33}3435bool garch_variance_path(std::span<const f64> r, std::span<const f64> th, GarchModel m, GarchDist d, std::vector<f64>& s2) {36 const usize n = r.size(); s2.assign(n, std::numeric_limits<f64>::quiet_NaN());37 if (n < 2) return false;38 const f64 mu = th[0], omega = th[1], alpha = th[2];39 const f64 gamma = m != GarchModel::garch ? th[3] : 0.0;40 const f64 beta = th[m != GarchModel::garch ? 4 : 3];41 // initial variance: sample variance of ε (ddof=0), as arch's default for the first observation uses backcast... we use42 // the plain sample variance — documented in the registry (tests pin this convention against arch with backcast given).43 f64 m1 = 0; for (f64 v : r) m1 += (v - mu); m1 /= static_cast<f64>(n);44 f64 var0 = 0; for (f64 v : r) var0 += (v - mu - m1) * (v - mu - m1); var0 /= static_cast<f64>(n);45 if (!(var0 > 0)) return false;46 if (m == GarchModel::egarch) {47 const f64 Eabs = (d == GarchDist::student_t) ? [&] { const f64 nu = th.back(); return 2.0 * std::sqrt(nu - 2.0) * std::exp(std::lgamma((nu + 1) / 2) - std::lgamma(nu / 2)) / ((nu - 1) * std::sqrt(std::numbers::pi)); }() : std::sqrt(2.0 / std::numbers::pi);48 f64 ls = std::log(var0);49 s2[0] = var0;50 for (usize t = 1; t < n; ++t) {51 const f64 e = r[t - 1] - mu, z = e / std::sqrt(s2[t - 1]);52 ls = omega + alpha * (std::abs(z) - Eabs) + gamma * z + beta * std::log(s2[t - 1]);53 s2[t] = std::exp(ls);54 if (!std::isfinite(s2[t]) || s2[t] <= 0) return false;55 }56 return true;57 }58 s2[0] = var0;59 for (usize t = 1; t < n; ++t) {60 const f64 e = r[t - 1] - mu; const f64 e2 = e * e;61 s2[t] = omega + (alpha + (e < 0 ? gamma : 0.0)) * e2 + beta * s2[t - 1];62 if (!std::isfinite(s2[t]) || s2[t] <= 0) return false;63 }64 return true;65}6667f64 garch_nll(std::span<const f64> r, std::span<const f64> th, GarchModel m, GarchDist d) {68 if (!garch_params_valid(th, m, d)) return std::numeric_limits<f64>::infinity();69 std::vector<f64> s2;70 if (!garch_variance_path(r, th, m, d, s2)) return std::numeric_limits<f64>::infinity();71 const f64 mu = th[0]; const usize n = r.size();72 f64 nll = 0, c = 0;73 auto add = [&](f64 v) { const f64 t = nll + v; if (std::abs(nll) >= std::abs(v)) c += (nll - t) + v; else c += (v - t) + nll; nll = t; };74 if (d == GarchDist::gaussian) {75 const f64 k = 0.5 * std::log(2.0 * std::numbers::pi);76 for (usize t = 0; t < n; ++t) { const f64 e = r[t] - mu; add(k + 0.5 * std::log(s2[t]) + 0.5 * e * e / s2[t]); }77 } else {78 const f64 nu = th.back();79 const f64 k = std::lgamma((nu + 1) / 2) - std::lgamma(nu / 2) - 0.5 * std::log(std::numbers::pi * (nu - 2));80 for (usize t = 0; t < n; ++t) { const f64 e = r[t] - mu; add(-(k - 0.5 * std::log(s2[t]) - 0.5 * (nu + 1) * std::log1p(e * e / ((nu - 2) * s2[t])))); }81 }82 return nll + c;83}8485f64 garch_nll(std::span<const f64> r, std::span<const f64> th, GarchModel m, f64 nu_fixed) {86 if (nu_fixed > 2.0) { std::vector<f64> t(th.begin(), th.end()); t.push_back(nu_fixed); return garch_nll(r, t, m, GarchDist::student_t); }87 return garch_nll(r, th, m, GarchDist::gaussian);88}8990// Transforms: u = [μ, log ω, logit-type for (α, β[, γ]), log(ν − 2)]91std::vector<f64> garch_to_unconstrained(std::span<const f64> th, GarchModel m, GarchDist d) {92 std::vector<f64> u(th.begin(), th.end());93 u[1] = m == GarchModel::egarch ? th[1] : std::log(th[1]);94 const bool g = m != GarchModel::garch; const usize ib = g ? 4 : 3;95 if (m == GarchModel::egarch) { u[2] = th[2]; u[3] = th[3]; u[4] = std::atanh(std::clamp(th[4], -0.999999, 0.999999)); }96 else {97 // persistence p = α + β (+γ/2) ∈ (0,1): u_p = logit(p); share s = α/p ∈ (0,1): u_s = logit(s); γ free for GJR via α+γ ≥ 0 → γ = exp(u_g) − α98 const f64 gamma = g ? th[3] : 0.0; const f64 p = th[2] + th[ib] + 0.5 * gamma; const f64 s = p > 0 ? th[2] / p : 0.5;99 u[2] = std::log(std::clamp(p, 1e-12, 1 - 1e-12) / (1 - std::clamp(p, 1e-12, 1 - 1e-12)));100 u[ib] = std::log(std::clamp(s, 1e-12, 1 - 1e-12) / (1 - std::clamp(s, 1e-12, 1 - 1e-12)));101 if (g) u[3] = std::log(std::max(th[2] + gamma, 1e-12));102 }103 if (d == GarchDist::student_t) { const f64 s = std::clamp((th.back() - kNuMin) / (kNuMax - kNuMin), 1e-9, 1 - 1e-9); u.back() = std::log(s / (1 - s)); }104 return u;105}106107std::vector<f64> garch_from_unconstrained(std::span<const f64> u, GarchModel m, GarchDist d) {108 std::vector<f64> th(u.begin(), u.end());109 th[1] = m == GarchModel::egarch ? u[1] : std::exp(u[1]);110 const bool g = m != GarchModel::garch;111 auto sig = [](f64 x) { return 1.0 / (1.0 + std::exp(-x)); };112 if (m == GarchModel::egarch) { th[2] = u[2]; th[3] = u[3]; th[4] = std::tanh(u[4]); }113 else if (!g) { const f64 p = sig(u[2]) * (1 - 1e-9), s = sig(u[3]); th[2] = s * p; th[3] = (1 - s) * p; }114 else {115 // GJR: p = α + β + γ/2, s = α/p → α = s p; ag = α + γ = exp(u_g) → γ = ag − α; β = p − α − γ/2 (must be ≥ 0: clamp via s)116 const f64 p = sig(u[2]) * (1 - 1e-9), s = sig(u[4]); const f64 alpha = s * p; const f64 ag = std::exp(u[3]); const f64 gamma = ag - alpha;117 f64 beta = p - alpha - 0.5 * gamma; if (beta < 0) { beta = 0; }118 th[2] = alpha; th[3] = gamma; th[4] = beta;119 }120 if (d == GarchDist::student_t) th.back() = kNuMin + (kNuMax - kNuMin) * sig(u.back());121 return th;122}123124} // namespace raid::models125