"""HAR-RV at the intraday scale (Corsi 2009 adapted): log RV over the next H minutes regressed on log trailing RV over
the last 30 minutes, 390 minutes (one day) and 5 days (1,950 minutes), plus time-of-day seasonal log variance.
Fitted by OLS on training origins; predictive distribution Student-t with zero mean, scale from exp(fitted) and a
degrees-of-freedom estimated from standardised training residuals (method of moments, clipped to [2.5, 50])."""
from __future__ import annotations

import numpy as np

from ..data.windows import SymbolGrid, N_DAY
from .naive import _hist_returns


def har_features(g: SymbolGrid, s: int, m: int, seas_prof: np.ndarray | None = None) -> np.ndarray:
    r, sl = _hist_returns(g, s, m, 6)
    if len(r) < 60:
        return None
    idx_today = np.flatnonzero(np.arange(len(r)) >= len(r) - m) if m > 0 else np.zeros(0, dtype=int)
    rv30 = np.mean(r[-30:] ** 2)
    rv390 = np.mean(r[-390:] ** 2)
    rv5d = np.mean(r ** 2)
    f = [np.log(rv30 + 1e-12), np.log(rv390 + 1e-12), np.log(rv5d + 1e-12)]
    return np.array(f)


class HARRV:
    name = "har_rv"

    def __init__(self):
        self.beta: dict[int, np.ndarray] = {}
        self.nu: dict[int, float] = {}
        self.sd: dict[int, float] = {}

    @staticmethod
    def _design(X: np.ndarray, m: np.ndarray, H: int) -> np.ndarray:
        tod = (m + 1 + H / 2) / 389.0
        return np.column_stack([np.ones(len(X)), X, tod, tod ** 2, np.log(H) * np.ones(len(X))])

    def fit(self, X: np.ndarray, m: np.ndarray, rv: np.ndarray, cum: np.ndarray, H: int) -> None:
        A = self._design(X, m, H)
        y = np.log(rv + 1e-12)
        beta, *_ = np.linalg.lstsq(A, y, rcond=None)
        self.beta[H] = beta
        pred_var = np.exp(A @ beta)
        z = cum / np.sqrt(pred_var)
        k = float(np.mean(z ** 4) / np.mean(z ** 2) ** 2)   # kurtosis -> nu = 4 + 6/(k-3)
        nu = 4 + 6 / max(k - 3, 0.12)
        self.nu[H] = float(np.clip(nu, 2.5, 50))
        self.sd[H] = float(np.std(y - A @ beta))

    def predict_batch(self, X: np.ndarray, m: np.ndarray, H: int) -> dict:
        A = self._design(X, m, H)
        var = np.exp(A @ self.beta[H] + 0.5 * self.sd[H] ** 2)   # log-normal mean correction
        nu = self.nu[H]
        return {"mu": np.zeros(len(X)), "scale": np.sqrt(var * (nu - 2) / nu), "nu": np.full(len(X), nu), "var": var}