1#include "raid/estimators/realized.hpp"2#include "raid/core/error.hpp"3#include "raid/statistics/reductions.hpp"4#include <cmath>5#include <numbers>67namespace raid::est {89f64 mu_p(f64 p) { return std::pow(2.0, p / 2.0) * std::tgamma((p + 1.0) / 2.0) / std::tgamma(0.5); }1011RealizedMeasures realized_measures_batch(const backend::Batch& b, const RealizedOptions& opt) {12 RealizedMeasures m; const usize S = b.n_series();13 std::string reason; auto be = backend::resolve(opt.plan, &reason);14 backend::KernelStats st1, st2, st3, st4; std::vector<f64> adj, s4; bool fused = false;15#ifdef RAID_ENABLE_METAL16 if (opt.fuse_kernels) fused = backend::realized_fused(*be, b, m.rv, adj, m.rsv_neg, m.rsv_pos, s4, &st1);17#endif18 if (!fused) {19 m.rv = be->sum_squares(b, &st1);20 adj = be->sum_abs_adjacent(b, &st2);21 be->semivariances(b, m.rsv_neg, m.rsv_pos, &st3);22 s4 = be->sum_fourth(b, &st4);23 }24 m.stats = st1; m.stats.kernel = fused ? "realized_fused" : "realized_measures(4 kernels)"; m.stats.wall_seconds = st1.wall_seconds + st2.wall_seconds + st3.wall_seconds + st4.wall_seconds;25 m.stats.device_seconds = st1.device_seconds + st2.device_seconds + st3.device_seconds + st4.device_seconds; m.stats.note = reason + (m.stats.note.empty() ? "" : " | " + m.stats.note);26 m.bpv.resize(S); m.rq.resize(S); m.jump.resize(S); m.tripv.resize(S); m.n_returns.resize(S); m.daily_return.resize(S);27 const f64 mu43 = mu_p(4.0 / 3.0);28 for (usize s = 0; s < S; ++s) {29 const usize M = b.size(s); m.n_returns[s] = M;30 auto x = b.series(s);31 m.daily_return[s] = stats::sum_neumaier(x);32 const f64 fs = (opt.bpv_finite_sample && M > 1) ? static_cast<f64>(M) / static_cast<f64>(M - 1) : 1.0;33 m.bpv[s] = M > 1 ? (1.0 / (kMu1 * kMu1)) * fs * adj[s] : std::numeric_limits<f64>::quiet_NaN();34 m.rq[s] = M > 0 ? (static_cast<f64>(M) / 3.0) * s4[s] : std::numeric_limits<f64>::quiet_NaN();35 m.jump[s] = std::max(m.rv[s] - m.bpv[s], 0.0);36 // tripower (CPU, FP64 — small cost, not batched on GPU in Tier A)37 if (M > 2) { f64 acc = 0, c = 0; for (usize j = 2; j < M; ++j) { const f64 v = std::pow(std::abs(x[j]), 4.0 / 3.0) * std::pow(std::abs(x[j - 1]), 4.0 / 3.0) * std::pow(std::abs(x[j - 2]), 4.0 / 3.0), t = acc + v; if (std::abs(acc) >= v) c += (acc - t) + v; else c += (v - t) + acc; acc = t; }38 const f64 fs3 = opt.bpv_finite_sample ? static_cast<f64>(M) / static_cast<f64>(M - 2) : 1.0; m.tripv[s] = fs3 * (acc + c) / (mu43 * mu43 * mu43); }39 else m.tripv[s] = std::numeric_limits<f64>::quiet_NaN();40 if (M < 2) m.warnings.push_back("segment " + std::to_string(s) + ": fewer than 2 returns");41 if (opt.annualize) { m.rv[s] *= opt.periods_per_year; m.bpv[s] *= opt.periods_per_year; m.rsv_neg[s] *= opt.periods_per_year; m.rsv_pos[s] *= opt.periods_per_year; m.jump[s] *= opt.periods_per_year; m.tripv[s] *= opt.periods_per_year; m.rq[s] *= opt.periods_per_year * opt.periods_per_year; }42 }43 return m;44}4546RealizedMeasures realized_measures(const data::IntradayReturns& ir, const RealizedOptions& opt) {47 backend::Batch b{ir.r, ir.session_offsets};48 auto m = realized_measures_batch(b, opt);49 m.sessions = ir.session_dates;50 return m;51}5253RealizedCov realized_covariance(const Eigen::Ref<const Eigen::MatrixXd>& R) {54 RAID_REQUIRE(R.rows() >= 1 && R.cols() >= 1, insufficient_samples, "realized_covariance: empty");55 RealizedCov out; out.n_obs = static_cast<usize>(R.rows());56 out.cov = R.transpose() * R; // Σ_j r_j r_j' (no demeaning: realized covariance of returns)57 const Eigen::Index N = R.cols(); out.corr.resize(N, N);58 for (Eigen::Index a = 0; a < N; ++a) for (Eigen::Index c = 0; c < N; ++c) { const f64 d = std::sqrt(out.cov(a, a) * out.cov(c, c)); out.corr(a, c) = d > 0 ? out.cov(a, c) / d : std::numeric_limits<f64>::quiet_NaN(); }59 return out;60}6162} // namespace raid::est63