// VAR(p), generalized FEVD (Pesaran–Shin 1998, Koop–Pesaran–Potter 1996) and Diebold–Yilmaz (2012, 2014) connectedness.
//   y_t = c + A_1 y_{t−1} + … + A_p y_{t−p} + u_t,  Σ = E[u u'] (OLS equation by equation = multivariate LS)
//   MA representation: y_t = Σ_h Φ_h u_{t−h}, Φ_0 = I, Φ_h = Σ_{j=1}^{min(h,p)} A_j Φ_{h−j}
//   GFEVD (horizon H): θ_{ij}(H) = σ_jj^{-1} Σ_{h=0}^{H−1} (e_i' Φ_h Σ e_j)² / Σ_{h=0}^{H−1} e_i' Φ_h Σ Φ_h' e_i,
//   normalised row-wise: θ̃_{ij} = θ_{ij} / Σ_j θ_{ij}  (Diebold–Yilmaz 2012 eq. 3–4)
//   Connectedness: total C = 100 · Σ_{i≠j} θ̃_{ij} / N ; from-others C_i← = Σ_{j≠i} θ̃_{ij} ; to-others C_→j = Σ_{i≠j} θ̃_{ij} ; net = to − from.
#pragma once
#include <vector>
#include <Eigen/Dense>
#include "raid/core/types.hpp"
#include "raid/core/estimate.hpp"

namespace raid::models {

struct VarFit {
  int p{1}; usize n_obs{0}; int N{0};
  Eigen::MatrixXd coef;        // N × (1 + N p): [c, A_1, …, A_p]
  Eigen::MatrixXd sigma;       // residual covariance (ddof = N p + 1, as statsmodels VAR `sigma_u`)
  Eigen::MatrixXd sigma_mle;   // ddof = 0
  Eigen::MatrixXd residuals;   // n_obs × N
  f64 loglik{0}, aic{0}, bic{0}, hqic{0};
  bool stable{true};           // all companion eigenvalues inside the unit circle
  f64 max_root{0};
  f64 rcond{0};
  std::vector<std::string> warnings;
};
/// Y is T × N (rows = time). Fits by OLS with an intercept (as statsmodels `VAR(Y).fit(p, trend="c")`).
VarFit var_fit(const Eigen::Ref<const Eigen::MatrixXd>& Y, int p);
/// Information criteria for p = 1..pmax on the common sample (observations pmax+1..T), as statsmodels `select_order`.
struct VarOrder { std::vector<int> p; std::vector<f64> aic, bic, hqic; int aic_best, bic_best, hqic_best; };
VarOrder var_select_order(const Eigen::Ref<const Eigen::MatrixXd>& Y, int pmax);
/// MA coefficient matrices Φ_0..Φ_{H−1}.
std::vector<Eigen::MatrixXd> var_ma(const VarFit& f, int H);
/// Generalized FEVD at horizon H (N × N, rows normalised to one).
Eigen::MatrixXd gfevd(const VarFit& f, int H);
/// Orthogonalised (Cholesky) FEVD at horizon H, for comparison with statsmodels `fevd`.
Eigen::MatrixXd fevd_cholesky(const VarFit& f, int H);

struct Connectedness { Eigen::MatrixXd table; f64 total; Eigen::VectorXd from, to, net; int H; };
Connectedness connectedness(const Eigen::MatrixXd& gfevd_table);
Connectedness connectedness(const VarFit& f, int H);

/// Rolling connectedness: windows of `window` rows, step `step`; one VAR(p) per window (CPU, parallel over windows).
struct RollingConnectedness { std::vector<usize> end_index; std::vector<f64> total; Eigen::MatrixXd net; usize n_windows{0}; usize failures{0}; f64 seconds{0}; };
RollingConnectedness rolling_connectedness(const Eigen::Ref<const Eigen::MatrixXd>& Y, int p, int H, usize window, usize step = 1);

/// Ledoit–Wolf (2004) linear shrinkage toward the scaled identity: Σ* = (1−δ) S + δ μ I, with the analytical δ*.
struct ShrunkCovariance { Eigen::MatrixXd cov; f64 shrinkage; f64 mu; };
ShrunkCovariance ledoit_wolf(const Eigen::Ref<const Eigen::MatrixXd>& X, bool assume_centered = false);

}  // namespace raid::models