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;15 m.rv = be->sum_squares(b, &st1);16 auto adj = be->sum_abs_adjacent(b, &st2);17 be->semivariances(b, m.rsv_neg, m.rsv_pos, &st3);18 auto s4 = be->sum_fourth(b, &st4);19 m.stats = st1; m.stats.kernel = "realized_measures"; m.stats.wall_seconds = st1.wall_seconds + st2.wall_seconds + st3.wall_seconds + st4.wall_seconds;20 m.stats.device_seconds = st1.device_seconds + st2.device_seconds + st3.device_seconds + st4.device_seconds; m.stats.note = reason;21 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);22 const f64 mu43 = mu_p(4.0 / 3.0);23 for (usize s = 0; s < S; ++s) {24 const usize M = b.size(s); m.n_returns[s] = M;25 auto x = b.series(s);26 m.daily_return[s] = stats::sum_neumaier(x);27 const f64 fs = (opt.bpv_finite_sample && M > 1) ? static_cast<f64>(M) / static_cast<f64>(M - 1) : 1.0;28 m.bpv[s] = M > 1 ? (1.0 / (kMu1 * kMu1)) * fs * adj[s] : std::numeric_limits<f64>::quiet_NaN();29 m.rq[s] = M > 0 ? (static_cast<f64>(M) / 3.0) * s4[s] : std::numeric_limits<f64>::quiet_NaN();30 m.jump[s] = std::max(m.rv[s] - m.bpv[s], 0.0);31 // tripower (CPU, FP64 — small cost, not batched on GPU in Tier A)32 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; }33 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); }34 else m.tripv[s] = std::numeric_limits<f64>::quiet_NaN();35 if (M < 2) m.warnings.push_back("segment " + std::to_string(s) + ": fewer than 2 returns");36 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; }37 }38 return m;39}4041RealizedMeasures realized_measures(const data::IntradayReturns& ir, const RealizedOptions& opt) {42 backend::Batch b{ir.r, ir.session_offsets};43 auto m = realized_measures_batch(b, opt);44 m.sessions = ir.session_dates;45 return m;46}4748RealizedCov realized_covariance(const Eigen::Ref<const Eigen::MatrixXd>& R) {49 RAID_REQUIRE(R.rows() >= 1 && R.cols() >= 1, insufficient_samples, "realized_covariance: empty");50 RealizedCov out; out.n_obs = static_cast<usize>(R.rows());51 out.cov = R.transpose() * R; // Σ_j r_j r_j' (no demeaning: realized covariance of returns)52 const Eigen::Index N = R.cols(); out.corr.resize(N, N);53 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(); }54 return out;55}5657} // namespace raid::est58