// GARCH-family conditional variance recursions and (quasi-)log-likelihoods. FP64 reference; the Metal kernel
// implements the same recursion in FP32 (documented tolerance). Definitions (registry key garch.*):
//   GARCH(1,1):     σ²_t = ω + α ε²_{t−1} + β σ²_{t−1}
//   GJR-GARCH(1,1): σ²_t = ω + (α + γ 1{ε_{t−1}<0}) ε²_{t−1} + β σ²_{t−1}
//   EGARCH(1,1):    ln σ²_t = ω + α (|z_{t−1}| − E|z|) + γ z_{t−1} + β ln σ²_{t−1},  z = ε/σ
// ε_t = r_t − μ. Initialisation: σ²_0 = sample variance of ε (GARCH/GJR), ln σ²_0 = ln(sample variance) (EGARCH) —
// the same convention as the `arch` package default ("backcast" is NOT used; see registry for the rugarch mapping).
// Densities: Gaussian, or standardized Student-t with ν > 2 (unit variance), as in arch/rugarch.
#pragma once
#include <span>
#include <vector>
#include <string>
#include "raid/core/types.hpp"
#include "raid/core/estimate.hpp"
#include "raid/core/plan.hpp"

namespace raid::models {

enum class GarchModel : int { garch = 0, gjr = 1, egarch = 2 };
enum class GarchDist : int { gaussian = 0, student_t = 1 };
std::string to_string(GarchModel m);
std::string to_string(GarchDist d);

/// Parameter layout θ (natural parameterisation): [μ, ω, α, (γ), β, (ν)] — γ only for gjr/egarch, ν only for student_t.
usize garch_n_params(GarchModel m, GarchDist d);
std::vector<std::string> garch_param_names(GarchModel m, GarchDist d);

/// Conditional variance path σ²_t (t = 0..n−1) for returns r and parameters θ; writes to sigma2 (resized).
/// Returns false (and leaves NaNs) if the recursion produced a non-positive or non-finite variance.
bool garch_variance_path(std::span<const f64> r, std::span<const f64> theta, GarchModel m, GarchDist d, std::vector<f64>& sigma2);

/// Negative average? NO: negative *total* log-likelihood Σ_t −log f(r_t | σ²_t). Returns +inf for invalid parameters
/// (ω ≤ 0, α < 0, β < 0, α+β ≥ 1 for GARCH; α+β+γ/2 ≥ 1 for GJR; |β| ≥ 1 for EGARCH; ν ≤ 2).
/// `nu_fixed`: if > 2 and the model is Gaussian-coded, a Student-t density with that ν is used (batched grids).
f64 garch_nll(std::span<const f64> r, std::span<const f64> theta, GarchModel m, f64 nu_fixed = 0.0);
f64 garch_nll(std::span<const f64> r, std::span<const f64> theta, GarchModel m, GarchDist d);

/// Stationarity / positivity constraints check.
bool garch_params_valid(std::span<const f64> theta, GarchModel m, GarchDist d);

/// Unconstrained ↔ constrained transforms used by the optimiser (ω = exp, α,β via logistic simplex, ν = 2 + exp).
std::vector<f64> garch_to_unconstrained(std::span<const f64> theta, GarchModel m, GarchDist d);
std::vector<f64> garch_from_unconstrained(std::span<const f64> u, GarchModel m, GarchDist d);

struct GarchFitOptions {
  GarchModel model{GarchModel::garch};
  GarchDist dist{GarchDist::gaussian};
  int max_iter{500};
  f64 tol{1e-8};              // relative change in NLL
  std::string cov{"qmle"};    // "qmle" (Bollerslev–Wooldridge sandwich) | "hessian" | "opg"
  int forecast_horizon{0};    // h ≥ 1 → multi-step variance forecasts in extra["forecast_variance"]
  bool scale_returns{true};   // internally scale returns to unit variance (as arch does, 'rescale')
  ExecutionPlan plan{};
};

/// Fit by (quasi-)maximum likelihood on the CPU (L-BFGS-style quasi-Newton on the unconstrained parameters with
/// numerical gradients validated against finite differences in tests). Never reports a non-converged fit as ok().
Estimate garch_fit(std::span<const f64> r, const GarchFitOptions& opt = {});
/// Multi-step-ahead conditional variance forecasts from the last state (h = 1..H). For EGARCH, h > 1 uses simulation-free
/// log-normal approximation only when `exact` is false; by default h > 1 forecasts for EGARCH are produced by
/// Monte Carlo with the given seed (reported in extra). GARCH/GJR forecasts are analytical.
std::vector<f64> garch_forecast(const Estimate& fit, std::span<const f64> r, int horizon, u64 seed = 20261102, int n_sim = 20000);

}  // namespace raid::models