// Robust covariance estimators for linear regression (CPU FP64 reference).
// Conventions follow statsmodels 0.15 (HC0–HC3, HAC Bartlett/Newey–West with use_correction=False) and
// Cameron–Gelbach–Miller one-way cluster with the Stata-style correction G/(G−1)·(N−1)/(N−K).
#pragma once
#include <vector>
#include <span>
#include <string>
#include <Eigen/Dense>
#include "raid/core/types.hpp"

namespace raid::stats {

enum class CovType { classical, hc0, hc1, hc2, hc3, hac_bartlett, cluster };
std::string to_string(CovType c);
CovType covtype_from_string(std::string_view s);

struct CovOptions {
  CovType type{CovType::classical};
  int hac_lags{-1};                 // −1 → floor(4 (n/100)^{2/9}) (Newey–West 1994 rule of thumb used by statsmodels users; documented)
  bool hac_small_sample{false};     // n/(n−k) factor (statsmodels use_correction)
  std::span<const i64> clusters{};  // cluster id per observation (CovType::cluster)
  bool cluster_small_sample{true};  // G/(G−1)·(N−1)/(N−K)
};

/// Covariance of OLS/WLS coefficients given X (n×k), residuals e (n), and (X'X)^{-1} computed stably by the caller.
/// `leverage` (h_ii) is required for HC2/HC3 and may be empty otherwise.
Eigen::MatrixXd robust_covariance(const Eigen::Ref<const Eigen::MatrixXd>& X, const Eigen::Ref<const Eigen::VectorXd>& resid,
                                  const Eigen::Ref<const Eigen::MatrixXd>& xtx_inv, const CovOptions& opt,
                                  const Eigen::VectorXd& leverage = Eigen::VectorXd());
int default_hac_lags(usize n);

}  // namespace raid::stats