1#include "raid/estimators/realized.hpp"2#include "raid/core/error.hpp"3#include "raid/statistics/reductions.hpp"4#include "raid/core/thread_pool.hpp"5#include <mutex>6#include <cmath>7#include <numbers>89namespace raid::est {1011f64 mu_p(f64 p) { return std::pow(2.0, p / 2.0) * std::tgamma((p + 1.0) / 2.0) / std::tgamma(0.5); }1213RealizedMeasures realized_measures_batch(const backend::Batch& b, const RealizedOptions& opt) {14 RealizedMeasures m; const usize S = b.n_series();15 std::string reason; auto be = backend::resolve(opt.plan, &reason);16 backend::KernelStats st1, st2, st3, st4; std::vector<f64> adj, s4; bool fused = false;17#ifdef RAID_ENABLE_METAL18 if (opt.fuse_kernels) fused = backend::realized_fused(*be, b, m.rv, adj, m.rsv_neg, m.rsv_pos, s4, &st1);19#endif20 if (!fused) {21 m.rv = be->sum_squares(b, &st1);22 adj = be->sum_abs_adjacent(b, &st2);23 be->semivariances(b, m.rsv_neg, m.rsv_pos, &st3);24 s4 = be->sum_fourth(b, &st4);25 }26 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;27 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);28 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);29 const f64 mu43 = mu_p(4.0 / 3.0);30 if (opt.tripower) {31 ThreadPool::global().parallel_for(S, [&](usize s0, usize s1) { std::vector<f64> p43; for (usize s = s0; s < s1; ++s) { const usize M = b.size(s); auto x = b.series(s);32 if (M <= 2) { m.tripv[s] = std::numeric_limits<f64>::quiet_NaN(); continue; }33 p43.resize(M); for (usize j = 0; j < M; ++j) { const f64 a = std::abs(x[j]); p43[j] = std::cbrt(a * a * a * a); } // |x|^{4/3} = cbrt(x⁴)34 f64 acc = 0, c = 0; for (usize j = 2; j < M; ++j) { const f64 v = p43[j] * p43[j - 1] * p43[j - 2], t = acc + v; if (std::abs(acc) >= v) c += (acc - t) + v; else c += (v - t) + acc; acc = t; }35 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); } });36 } else m.tripv.assign(S, std::numeric_limits<f64>::quiet_NaN());37 std::mutex warn_m;38 ThreadPool::global().parallel_for(S, [&](usize s0, usize s1) { for (usize s = s0; s < s1; ++s) {39 const usize M = b.size(s); m.n_returns[s] = M;40 auto x = b.series(s);41 m.daily_return[s] = stats::sum_neumaier(x);42 const f64 fs = (opt.bpv_finite_sample && M > 1) ? static_cast<f64>(M) / static_cast<f64>(M - 1) : 1.0;43 m.bpv[s] = M > 1 ? (1.0 / (kMu1 * kMu1)) * fs * adj[s] : std::numeric_limits<f64>::quiet_NaN();44 m.rq[s] = M > 0 ? (static_cast<f64>(M) / 3.0) * s4[s] : std::numeric_limits<f64>::quiet_NaN();45 m.jump[s] = std::max(m.rv[s] - m.bpv[s], 0.0);46 if (M < 2) { std::lock_guard<std::mutex> lk(warn_m); m.warnings.push_back("segment " + std::to_string(s) + ": fewer than 2 returns"); }47 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; }48 } });49 return m;50}5152RealizedMeasures realized_measures(const data::IntradayReturns& ir, const RealizedOptions& opt) {53 backend::Batch b{ir.r, ir.session_offsets};54 auto m = realized_measures_batch(b, opt);55 m.sessions = ir.session_dates;56 return m;57}5859RealizedCov realized_covariance(const Eigen::Ref<const Eigen::MatrixXd>& R) {60 RAID_REQUIRE(R.rows() >= 1 && R.cols() >= 1, insufficient_samples, "realized_covariance: empty");61 RealizedCov out; out.n_obs = static_cast<usize>(R.rows());62 out.cov = R.transpose() * R; // Σ_j r_j r_j' (no demeaning: realized covariance of returns)63 const Eigen::Index N = R.cols(); out.corr.resize(N, N);64 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(); }65 return out;66}6768} // namespace raid::est69