1#include "raid/models/var.hpp"2#include "raid/core/error.hpp"3#include "raid/core/thread_pool.hpp"4#include <Eigen/Eigenvalues>5#include <chrono>6#include <cmath>7#include <numbers>89namespace raid::models {1011namespace {12void build_design(const Eigen::Ref<const Eigen::MatrixXd>& Y, int p, usize start, Eigen::MatrixXd& X, Eigen::MatrixXd& Yt) {13 const Eigen::Index T = Y.rows(), N = Y.cols(); const Eigen::Index n = T - static_cast<Eigen::Index>(start);14 X.resize(n, 1 + N * p); Yt.resize(n, N);15 for (Eigen::Index i = 0; i < n; ++i) { const Eigen::Index t = static_cast<Eigen::Index>(start) + i; X(i, 0) = 1.0; for (int l = 1; l <= p; ++l) X.block(i, 1 + (l - 1) * N, 1, N) = Y.row(t - l); Yt.row(i) = Y.row(t); }16}17}1819VarFit var_fit(const Eigen::Ref<const Eigen::MatrixXd>& Y, int p) {20 const Eigen::Index T = Y.rows(), N = Y.cols(); RAID_REQUIRE(p >= 1 && T > p + N * p + 2, insufficient_samples, "var_fit: need T > p + N·p + 2");21 RAID_REQUIRE(Y.allFinite(), data_quality, "var_fit: non-finite values");22 VarFit f; f.p = p; f.N = static_cast<int>(N); Eigen::MatrixXd X, Yt; build_design(Y, p, static_cast<usize>(p), X, Yt); f.n_obs = static_cast<usize>(X.rows());23 Eigen::ColPivHouseholderQR<Eigen::MatrixXd> qr(X); qr.setThreshold(1e-12); const Eigen::Index k = X.cols();24 const Eigen::MatrixXd R = qr.matrixR().topLeftCorner(k, k).template triangularView<Eigen::Upper>(); f.rcond = std::abs(R(k - 1, k - 1)) / std::abs(R(0, 0));25 if (qr.rank() < k) f.warnings.push_back("VAR design rank deficient (rank " + std::to_string(qr.rank()) + " < " + std::to_string(k) + ")");26 Eigen::MatrixXd B = qr.solve(Yt); // k × N27 f.coef = B.transpose(); f.residuals = Yt - X * B; const f64 n = static_cast<f64>(f.n_obs);28 f.sigma_mle = f.residuals.transpose() * f.residuals / n; f.sigma = f.residuals.transpose() * f.residuals / (n - static_cast<f64>(k));29 // log-likelihood (Gaussian, MLE covariance), as statsmodels VARResults.llf30 const f64 ld = std::log(f.sigma_mle.determinant()); const f64 Nd = static_cast<f64>(N);31 f.loglik = -0.5 * n * (Nd * std::log(2 * std::numbers::pi) + ld + Nd);32 const f64 nparams = static_cast<f64>(k) * Nd; // statsmodels: ln|Σ_mle| + 2 k N / n etc. (uses sigma with nobs correction: see select_order); we report the mle-based criteria33 f.aic = ld + 2.0 * nparams / n; f.bic = ld + std::log(n) * nparams / n; f.hqic = ld + 2.0 * std::log(std::log(n)) * nparams / n;34 // stability: companion matrix eigenvalues35 Eigen::MatrixXd C = Eigen::MatrixXd::Zero(N * p, N * p); for (int l = 1; l <= p; ++l) C.block(0, (l - 1) * N, N, N) = f.coef.block(0, 1 + (l - 1) * N, N, N); if (p > 1) C.block(N, 0, N * (p - 1), N * (p - 1)) = Eigen::MatrixXd::Identity(N * (p - 1), N * (p - 1));36 Eigen::EigenSolver<Eigen::MatrixXd> es(C, false); f.max_root = es.eigenvalues().cwiseAbs().maxCoeff(); f.stable = f.max_root < 1.0; if (!f.stable) f.warnings.push_back("VAR not stable (max root " + std::to_string(f.max_root) + ")");37 return f;38}3940VarOrder var_select_order(const Eigen::Ref<const Eigen::MatrixXd>& Y, int pmax) {41 VarOrder o; const Eigen::Index N = Y.cols(); const f64 Nd = static_cast<f64>(N);42 for (int p = 1; p <= pmax; ++p) { Eigen::MatrixXd X, Yt; build_design(Y, p, static_cast<usize>(pmax), X, Yt); Eigen::ColPivHouseholderQR<Eigen::MatrixXd> qr(X); Eigen::MatrixXd B = qr.solve(Yt); Eigen::MatrixXd U = Yt - X * B; const f64 n = static_cast<f64>(X.rows()); const f64 ld = std::log((U.transpose() * U / n).determinant()); const f64 np = static_cast<f64>(X.cols()) * Nd;43 o.p.push_back(p); o.aic.push_back(ld + 2 * np / n); o.bic.push_back(ld + std::log(n) * np / n); o.hqic.push_back(ld + 2 * std::log(std::log(n)) * np / n); }44 auto best = [&](const std::vector<f64>& v) { int b = 0; for (usize i = 1; i < v.size(); ++i) if (v[i] < v[static_cast<usize>(b)]) b = static_cast<int>(i); return o.p[static_cast<usize>(b)]; };45 o.aic_best = best(o.aic); o.bic_best = best(o.bic); o.hqic_best = best(o.hqic); return o;46}4748std::vector<Eigen::MatrixXd> var_ma(const VarFit& f, int H) {49 const int N = f.N; std::vector<Eigen::MatrixXd> Phi(static_cast<usize>(H)); Phi[0] = Eigen::MatrixXd::Identity(N, N);50 for (int h = 1; h < H; ++h) { Phi[static_cast<usize>(h)] = Eigen::MatrixXd::Zero(N, N); for (int j = 1; j <= std::min(h, f.p); ++j) Phi[static_cast<usize>(h)] += f.coef.block(0, 1 + (j - 1) * N, N, N) * Phi[static_cast<usize>(h - j)]; }51 return Phi;52}5354Eigen::MatrixXd gfevd(const VarFit& f, int H) {55 RAID_REQUIRE(H >= 1, invalid_argument, "gfevd: H >= 1"); const int N = f.N; auto Phi = var_ma(f, H); const Eigen::MatrixXd& S = f.sigma;56 Eigen::MatrixXd num = Eigen::MatrixXd::Zero(N, N); Eigen::VectorXd den = Eigen::VectorXd::Zero(N);57 for (int h = 0; h < H; ++h) { const Eigen::MatrixXd PS = Phi[static_cast<usize>(h)] * S; for (int i = 0; i < N; ++i) { den[i] += (PS.row(i) * Phi[static_cast<usize>(h)].row(i).transpose())(0, 0); for (int j = 0; j < N; ++j) num(i, j) += PS(i, j) * PS(i, j) / S(j, j); } }58 Eigen::MatrixXd th(N, N); for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) th(i, j) = num(i, j) / den[i];59 for (int i = 0; i < N; ++i) th.row(i) /= th.row(i).sum();60 return th;61}6263Eigen::MatrixXd fevd_cholesky(const VarFit& f, int H) {64 const int N = f.N; auto Phi = var_ma(f, H); Eigen::LLT<Eigen::MatrixXd> llt(f.sigma); RAID_REQUIRE(llt.info() == Eigen::Success, numerical, "fevd: Σ not PD"); const Eigen::MatrixXd P = llt.matrixL();65 Eigen::MatrixXd num = Eigen::MatrixXd::Zero(N, N); Eigen::VectorXd den = Eigen::VectorXd::Zero(N);66 for (int h = 0; h < H; ++h) { const Eigen::MatrixXd Th = Phi[static_cast<usize>(h)] * P; for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) { num(i, j) += Th(i, j) * Th(i, j); den[i] += Th(i, j) * Th(i, j); } }67 Eigen::MatrixXd out(N, N); for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) out(i, j) = num(i, j) / den[i]; return out;68}6970Connectedness connectedness(const Eigen::MatrixXd& th) {71 const Eigen::Index N = th.rows(); Connectedness c; c.table = th; c.from.resize(N); c.to.resize(N); c.net.resize(N); f64 off = 0;72 for (Eigen::Index i = 0; i < N; ++i) { f64 fr = 0, to = 0; for (Eigen::Index j = 0; j < N; ++j) if (i != j) { fr += th(i, j); to += th(j, i); off += th(i, j); } c.from[i] = 100.0 * fr; c.to[i] = 100.0 * to; c.net[i] = c.to[i] - c.from[i]; }73 c.total = 100.0 * off / static_cast<f64>(N); c.H = 0; return c;74}75Connectedness connectedness(const VarFit& f, int H) { auto c = connectedness(gfevd(f, H)); c.H = H; return c; }7677RollingConnectedness rolling_connectedness(const Eigen::Ref<const Eigen::MatrixXd>& Y, int p, int H, usize window, usize step) {78 const auto t0 = std::chrono::steady_clock::now(); RollingConnectedness r; const usize T = static_cast<usize>(Y.rows()), N = static_cast<usize>(Y.cols());79 RAID_REQUIRE(window <= T && window > static_cast<usize>(p + static_cast<int>(N) * p + 10), insufficient_samples, "rolling_connectedness: window");80 const usize nw = (T - window) / step + 1; r.end_index.resize(nw); r.total.assign(nw, std::numeric_limits<f64>::quiet_NaN()); r.net.resize(static_cast<Eigen::Index>(N), static_cast<Eigen::Index>(nw)); r.net.setConstant(std::numeric_limits<f64>::quiet_NaN()); std::atomic<usize> fails{0};81 Eigen::MatrixXd Yc = Y;82 ThreadPool::global().parallel_for(nw, [&](usize a, usize b) { for (usize w = a; w < b; ++w) { const usize beg = w * step; r.end_index[w] = beg + window; try { auto f = var_fit(Yc.middleRows(static_cast<Eigen::Index>(beg), static_cast<Eigen::Index>(window)), p); auto c = connectedness(f, H); r.total[w] = c.total; r.net.col(static_cast<Eigen::Index>(w)) = c.net; } catch (const RaidError&) { ++fails; } } });83 r.n_windows = nw; r.failures = fails.load(); r.seconds = std::chrono::duration<f64>(std::chrono::steady_clock::now() - t0).count(); return r;84}8586ShrunkCovariance ledoit_wolf(const Eigen::Ref<const Eigen::MatrixXd>& X0, bool assume_centered) {87 // Ledoit & Wolf (2004) "A well-conditioned estimator for large-dimensional covariance matrices", as sklearn.covariance.ledoit_wolf88 const Eigen::Index n = X0.rows(), p = X0.cols(); RAID_REQUIRE(n >= 2, insufficient_samples, "ledoit_wolf: n >= 2");89 Eigen::MatrixXd X = X0; if (!assume_centered) X.rowwise() -= X0.colwise().mean();90 const Eigen::MatrixXd S = X.transpose() * X / static_cast<f64>(n); const f64 mu = S.trace() / static_cast<f64>(p);91 // δ² = ||S − μI||²_F / p ; β̄ = (1/n²) Σ_i ||x_i x_i' − S||²_F / p (sklearn formula); shrinkage = min(β̄/δ², 1)92 const f64 delta2 = (S - mu * Eigen::MatrixXd::Identity(p, p)).squaredNorm() / static_cast<f64>(p);93 f64 beta = 0; for (Eigen::Index i = 0; i < n; ++i) { const Eigen::MatrixXd xx = X.row(i).transpose() * X.row(i); beta += (xx - S).squaredNorm(); } beta /= static_cast<f64>(n) * static_cast<f64>(n) * static_cast<f64>(p);94 const f64 shrink = delta2 > 0 ? std::min(1.0, beta / delta2) : 0.0;95 ShrunkCovariance out; out.shrinkage = shrink; out.mu = mu; out.cov = (1 - shrink) * S + shrink * mu * Eigen::MatrixXd::Identity(p, p); return out;96}9798} // namespace raid::models99