// HAR family (registry keys har.*). Daily series RV_t (t = 1..T), trading-day lags, causal construction:
//   HAR-RV   (Corsi 2009):  y_{t+h} = β0 + βd RV_t + βw RV^{(5)}_t + βm RV^{(22)}_t + ε,  RV^{(k)}_t = (1/k) Σ_{j=0}^{k−1} RV_{t−j}
//   log-HAR:               same with log RV on both sides (forecast = exp(ŷ) × smearing optional, not applied by default)
//   HAR-X:                 + user covariates X_t (already aligned at t)
//   HARQ (Bollerslev–Patton–Quaedvlieg 2016): βd → (βd + βq √RQ_t): term βq · (√RQ_t − mean √RQ) · RV_t added (demeaning optional)
//   Target for horizon h: y_{t+h} = (1/h) Σ_{j=1}^{h} RV_{t+j} (average over the next h days, "overlapping" when h > 1)
//   or the single-day RV_{t+h} (`target = single`). Overlapping targets induce MA(h−1) errors → HAC with ≥ h−1 lags.
#pragma once
#include <vector>
#include <span>
#include <string>
#include <Eigen/Dense>
#include "raid/core/types.hpp"
#include "raid/core/estimate.hpp"
#include "raid/core/plan.hpp"
#include "raid/estimators/ols.hpp"
#include "raid/backend/backend.hpp"

namespace raid::est {

enum class HarTarget { average, single };

struct HarSpec {
  std::vector<usize> lags{1, 5, 22};   // daily, weekly, monthly windows
  int horizon{1};
  HarTarget target{HarTarget::average};
  bool log{false};                      // log-HAR
  bool harq{false};                     // add √RQ_t · RV_t
  bool harq_demean_sqrt_rq{true};       // use (√RQ_t − √(mean RQ)) · RV_t exactly as highfrequency::HARmodel(type="HARQ"); the mean is over the sample passed to the fit (estimation window), never the future
  stats::CovOptions cov{stats::CovType::hac_bartlett, -1, false, {}, true};  // hac_lags −1 → max(h−1, ⌊4(n/100)^{2/9}⌋)
  std::string model_key() const { return std::string(log ? "har.logrv" : "har.rv") + (harq ? ".harq" : "") + ".h" + std::to_string(horizon); }
};

/// Builds the aligned design (y, X) from a daily RV series (and optional RQ and exogenous columns X_t, aligned at t).
struct HarDesign { Eigen::MatrixXd X; Eigen::VectorXd y; std::vector<std::string> names; usize first_t{0}; usize n{0}; f64 sqrt_rq_mean{0}; };
HarDesign har_design(std::span<const f64> rv, const HarSpec& spec, std::span<const f64> rq = {}, const Eigen::MatrixXd& exog = Eigen::MatrixXd(), const std::vector<std::string>& exog_names = {});

/// Fit one HAR on one series (CPU QR + robust covariance). Estimate.extra contains fitted/residual summaries.
Estimate har_fit(std::span<const f64> rv, const HarSpec& spec, std::span<const f64> rq = {}, const Eigen::MatrixXd& exog = Eigen::MatrixXd(), const std::vector<std::string>& exog_names = {});

/// Forecast ŷ for the last available t (origin = last observation) from a fitted HAR; for log models returns exp(ŷ)
/// (naive back-transform; the smearing factor is reported in extra["smearing"] and not applied).
f64 har_forecast(const Estimate& fit, std::span<const f64> rv, const HarSpec& spec, std::span<const f64> rq = {}, const Eigen::VectorXd& exog_last = Eigen::VectorXd());

/// Rolling (or expanding when window = 0) re-estimation with one-origin-ahead forecasts, strictly causal:
/// the model fitted on observations ≤ t predicts y_{t+h}. Returns forecasts aligned with the target index.
struct HarRolling { std::vector<usize> origin; std::vector<f64> forecast, actual; Eigen::MatrixXd coefficients; usize n_fits{0}; f64 seconds{0}; };
HarRolling har_rolling(std::span<const f64> rv, const HarSpec& spec, usize window, usize step = 1, std::span<const f64> rq = {});

/// Batched HAR over a panel of daily RV series (CSR). Sufficient statistics are accumulated by the backend (Metal when
/// available: regression_stats) and solved in FP64 on the CPU; inference is classical (HAC requires residuals: use
/// har_fit per series when robust SE are needed). Returns one Estimate per series.
struct HarBatchResult { std::vector<Estimate> fits; backend::KernelStats stats; };
HarBatchResult har_batch(std::span<const f64> rv_concat, std::span<const i64> offsets, const HarSpec& spec, const ExecutionPlan& plan = {});

/// Forecast losses for volatility: MSE and QLIKE (Patton 2011): QLIKE = y/ŷ − log(y/ŷ) − 1.
struct ForecastLoss { f64 mse, qlike; usize n; };
ForecastLoss forecast_loss(std::span<const f64> actual, std::span<const f64> forecast);

}  // namespace raid::est