"""Outils économétriques : hédonique semi-log, indices (dummies temporelles, ventes répétées BMN / Case–Shiller),
statistiques de ratio IAAO (COD, PRD, PRB), autocorrélation spatiale (I de Moran), matrices de poids kNN/distance."""
from __future__ import annotations
import numpy as np

# ----------------------------------------------------------------------------- MCO

def mco(X: np.ndarray, y: np.ndarray, robuste: bool = True) -> dict:
    """Moindres carrés ordinaires avec erreurs-types classiques et HC1 (White). X doit contenir la constante."""
    X = np.asarray(X, float); y = np.asarray(y, float)
    n, k = X.shape
    XtX_inv = np.linalg.pinv(X.T @ X)
    beta = XtX_inv @ X.T @ y
    res = y - X @ beta
    s2 = res @ res / (n - k)
    V_cl = s2 * XtX_inv
    meat = (X * (res[:, None] ** 2)).T @ X
    V_hc1 = XtX_inv @ meat @ XtX_inv * n / (n - k)
    tss = ((y - y.mean()) ** 2).sum()
    r2 = 1 - (res @ res) / tss
    return {"beta": beta, "se": np.sqrt(np.diag(V_cl)), "se_hc1": np.sqrt(np.diag(V_hc1)), "residus": res,
            "sigma2": s2, "r2": r2, "r2_ajuste": 1 - (1 - r2) * (n - 1) / (n - k), "n": n, "k": k, "V_hc1": V_hc1}

def indice_dummies_temps(beta_temps: np.ndarray, base: float = 100.0) -> np.ndarray:
    """Indice hédonique à partir des coefficients des indicatrices de période (période 0 = référence, coefficient 0)."""
    b = np.concatenate([[0.0], np.asarray(beta_temps, float)])
    return base * np.exp(b)

# ----------------------------------------------------------------------------- ventes répétées

def ventes_repetees_bmn(prix_1, prix_2, t_1, t_2, n_periodes: int, base: float = 100.0) -> dict:
    """Estimateur de Bailey–Muth–Nourse (1963) : ln(P2/P1) = Σ_t D_t β_t + ε, D = −1 à t1, +1 à t2, période 0 omise."""
    p1 = np.asarray(prix_1, float); p2 = np.asarray(prix_2, float)
    t1 = np.asarray(t_1, int); t2 = np.asarray(t_2, int)
    y = np.log(p2 / p1)
    D = np.zeros((len(y), n_periodes))
    for i in range(len(y)):
        D[i, t1[i]] -= 1; D[i, t2[i]] += 1
    D = D[:, 1:]
    beta, *_ = np.linalg.lstsq(D, y, rcond=None)
    res = y - D @ beta
    idx = base * np.exp(np.concatenate([[0.0], beta]))
    return {"beta": beta, "indice": idx, "residus": res, "D": D, "y": y}

def ventes_repetees_case_shiller(prix_1, prix_2, t_1, t_2, n_periodes: int, base: float = 100.0) -> dict:
    """Case–Shiller (1987, 1989) pondéré en trois étapes : (1) BMN, (2) régression de ε² sur l'intervalle,
    (3) moindres carrés pondérés par 1/√(variance ajustée)."""
    e1 = ventes_repetees_bmn(prix_1, prix_2, t_1, t_2, n_periodes, base)
    gap = (np.asarray(t_2) - np.asarray(t_1)).astype(float)
    Z = np.column_stack([np.ones_like(gap), gap])
    a, *_ = np.linalg.lstsq(Z, e1["residus"] ** 2, rcond=None)
    var_hat = np.maximum(Z @ a, 1e-8)
    w = 1 / np.sqrt(var_hat)
    beta, *_ = np.linalg.lstsq(e1["D"] * w[:, None], e1["y"] * w, rcond=None)
    idx = base * np.exp(np.concatenate([[0.0], beta]))
    return {"beta": beta, "indice": idx, "coef_variance": a, "poids": w}

# ----------------------------------------------------------------------------- ratios IAAO

def statistiques_ratio(valeurs_estimees, prix) -> dict:
    """Statistiques d'étude de ratios (IAAO, Standard on Ratio Studies) :
    ratio médian, COD = 100 × moyenne|r − médiane| / médiane, PRD = moyenne(r) / (Σ V / Σ P),
    PRB = pente de (r − médiane)/médiane sur log2 de la valeur proxy ((V/médiane + P)/2)."""
    V = np.asarray(valeurs_estimees, float); P = np.asarray(prix, float)
    r = V / P
    med = np.median(r)
    cod = 100.0 * np.mean(np.abs(r - med)) / med
    prd = np.mean(r) / (V.sum() / P.sum())
    proxy = 0.5 * (V / med + P)
    x = np.log2(proxy)
    yv = (r - med) / med
    X = np.column_stack([np.ones_like(x), x])
    b, *_ = np.linalg.lstsq(X, yv, rcond=None)
    res = yv - X @ b
    se = np.sqrt((res @ res) / (len(x) - 2) / ((x - x.mean()) ** 2).sum())
    return {"n": len(r), "ratio_median": med, "ratio_moyen": r.mean(), "COD": cod, "PRD": prd, "PRB": b[1], "PRB_se": se,
            "ratio_pondere": V.sum() / P.sum()}

# ----------------------------------------------------------------------------- spatial

def matrice_poids_knn(coords: np.ndarray, k: int = 8, standardiser: bool = True) -> np.ndarray:
    """Matrice W (n×n) des k plus proches voisins (distance euclidienne), standardisée en ligne."""
    c = np.asarray(coords, float)
    d = np.sqrt(((c[:, None, :] - c[None, :, :]) ** 2).sum(-1))
    np.fill_diagonal(d, np.inf)
    W = np.zeros_like(d)
    idx = np.argsort(d, axis=1)[:, :k]
    for i in range(len(c)):
        W[i, idx[i]] = 1.0
    if standardiser:
        W = W / W.sum(1, keepdims=True)
    return W

def matrice_poids_distance_inverse(coords: np.ndarray, portee: float, puissance: float = 1.0, standardiser: bool = True) -> np.ndarray:
    c = np.asarray(coords, float)
    d = np.sqrt(((c[:, None, :] - c[None, :, :]) ** 2).sum(-1))
    with np.errstate(divide="ignore"):
        W = np.where((d > 0) & (d <= portee), 1.0 / d ** puissance, 0.0)
    if standardiser:
        s = W.sum(1, keepdims=True); s[s == 0] = 1
        W = W / s
    return W

def moran_i(x: np.ndarray, W: np.ndarray) -> dict:
    """I de Moran et test sous hypothèse de normalité : I = (n/S0) · (z'Wz)/(z'z)."""
    z = np.asarray(x, float); z = z - z.mean()
    n = len(z); S0 = W.sum()
    I = (n / S0) * (z @ W @ z) / (z @ z)
    EI = -1.0 / (n - 1)
    S1 = 0.5 * ((W + W.T) ** 2).sum()
    S2 = ((W.sum(0) + W.sum(1)) ** 2).sum()
    VI = (n ** 2 * S1 - n * S2 + 3 * S0 ** 2) / ((n ** 2 - 1) * S0 ** 2) - EI ** 2
    zI = (I - EI) / np.sqrt(VI)
    from math import erf, sqrt
    p = 2 * (1 - 0.5 * (1 + erf(abs(zI) / sqrt(2))))
    return {"I": I, "E_I": EI, "Var_I": VI, "z": zI, "p": p}

def decalage_spatial(W: np.ndarray, x: np.ndarray) -> np.ndarray:
    return W @ np.asarray(x, float)

def sar_lag_mv(y: np.ndarray, X: np.ndarray, W: np.ndarray, grille_rho=None) -> dict:
    """Modèle autorégressif spatial y = ρWy + Xβ + ε estimé par maximum de vraisemblance concentrée (recherche sur ρ)."""
    y = np.asarray(y, float); X = np.asarray(X, float); n = len(y)
    Wy = W @ y
    b0, *_ = np.linalg.lstsq(X, y, rcond=None); e0 = y - X @ b0
    bL, *_ = np.linalg.lstsq(X, Wy, rcond=None); eL = Wy - X @ bL
    eig = np.linalg.eigvals(W).real
    rhos = grille_rho if grille_rho is not None else np.linspace(-0.9, 0.99, 380)
    best = None
    for rho in rhos:
        e = e0 - rho * eL
        s2 = e @ e / n
        ll = np.log(np.abs(1 - rho * eig)).sum() - n / 2 * np.log(2 * np.pi * s2) - n / 2
        if best is None or ll > best[0]:
            best = (ll, rho, s2)
    ll, rho, s2 = best
    beta = b0 - rho * bL
    return {"rho": rho, "beta": beta, "sigma2": s2, "loglik": ll}