spb/traite-evaluation-qc

1"""Outils économétriques : hédonique semi-log, indices (dummies temporelles, ventes répétées BMN / Case–Shiller),2statistiques de ratio IAAO (COD, PRD, PRB), autocorrélation spatiale (I de Moran), matrices de poids kNN/distance."""3from __future__ import annotations4import numpy as np56# ----------------------------------------------------------------------------- MCO78def mco(X: np.ndarray, y: np.ndarray, robuste: bool = True) -> dict:9    """Moindres carrés ordinaires avec erreurs-types classiques et HC1 (White). X doit contenir la constante."""10    X = np.asarray(X, float); y = np.asarray(y, float)11    n, k = X.shape12    XtX_inv = np.linalg.pinv(X.T @ X)13    beta = XtX_inv @ X.T @ y14    res = y - X @ beta15    s2 = res @ res / (n - k)16    V_cl = s2 * XtX_inv17    meat = (X * (res[:, None] ** 2)).T @ X18    V_hc1 = XtX_inv @ meat @ XtX_inv * n / (n - k)19    tss = ((y - y.mean()) ** 2).sum()20    r2 = 1 - (res @ res) / tss21    return {"beta": beta, "se": np.sqrt(np.diag(V_cl)), "se_hc1": np.sqrt(np.diag(V_hc1)), "residus": res,22            "sigma2": s2, "r2": r2, "r2_ajuste": 1 - (1 - r2) * (n - 1) / (n - k), "n": n, "k": k, "V_hc1": V_hc1}2324def indice_dummies_temps(beta_temps: np.ndarray, base: float = 100.0) -> np.ndarray:25    """Indice hédonique à partir des coefficients des indicatrices de période (période 0 = référence, coefficient 0)."""26    b = np.concatenate([[0.0], np.asarray(beta_temps, float)])27    return base * np.exp(b)2829# ----------------------------------------------------------------------------- ventes répétées3031def ventes_repetees_bmn(prix_1, prix_2, t_1, t_2, n_periodes: int, base: float = 100.0) -> dict:32    """Estimateur de Bailey–Muth–Nourse (1963) : ln(P2/P1) = Σ_t D_t β_t + ε, D = −1 à t1, +1 à t2, période 0 omise."""33    p1 = np.asarray(prix_1, float); p2 = np.asarray(prix_2, float)34    t1 = np.asarray(t_1, int); t2 = np.asarray(t_2, int)35    y = np.log(p2 / p1)36    D = np.zeros((len(y), n_periodes))37    for i in range(len(y)):38        D[i, t1[i]] -= 1; D[i, t2[i]] += 139    D = D[:, 1:]40    beta, *_ = np.linalg.lstsq(D, y, rcond=None)41    res = y - D @ beta42    idx = base * np.exp(np.concatenate([[0.0], beta]))43    return {"beta": beta, "indice": idx, "residus": res, "D": D, "y": y}4445def ventes_repetees_case_shiller(prix_1, prix_2, t_1, t_2, n_periodes: int, base: float = 100.0) -> dict:46    """Case–Shiller (1987, 1989) pondéré en trois étapes : (1) BMN, (2) régression de ε² sur l'intervalle,47    (3) moindres carrés pondérés par 1/√(variance ajustée)."""48    e1 = ventes_repetees_bmn(prix_1, prix_2, t_1, t_2, n_periodes, base)49    gap = (np.asarray(t_2) - np.asarray(t_1)).astype(float)50    Z = np.column_stack([np.ones_like(gap), gap])51    a, *_ = np.linalg.lstsq(Z, e1["residus"] ** 2, rcond=None)52    var_hat = np.maximum(Z @ a, 1e-8)53    w = 1 / np.sqrt(var_hat)54    beta, *_ = np.linalg.lstsq(e1["D"] * w[:, None], e1["y"] * w, rcond=None)55    idx = base * np.exp(np.concatenate([[0.0], beta]))56    return {"beta": beta, "indice": idx, "coef_variance": a, "poids": w}5758# ----------------------------------------------------------------------------- ratios IAAO5960def statistiques_ratio(valeurs_estimees, prix) -> dict:61    """Statistiques d'étude de ratios (IAAO, Standard on Ratio Studies) :62    ratio médian, COD = 100 × moyenne|r − médiane| / médiane, PRD = moyenne(r) / (Σ V / Σ P),63    PRB = pente de (r − médiane)/médiane sur log2 de la valeur proxy ((V/médiane + P)/2)."""64    V = np.asarray(valeurs_estimees, float); P = np.asarray(prix, float)65    r = V / P66    med = np.median(r)67    cod = 100.0 * np.mean(np.abs(r - med)) / med68    prd = np.mean(r) / (V.sum() / P.sum())69    proxy = 0.5 * (V / med + P)70    x = np.log2(proxy)71    yv = (r - med) / med72    X = np.column_stack([np.ones_like(x), x])73    b, *_ = np.linalg.lstsq(X, yv, rcond=None)74    res = yv - X @ b75    se = np.sqrt((res @ res) / (len(x) - 2) / ((x - x.mean()) ** 2).sum())76    return {"n": len(r), "ratio_median": med, "ratio_moyen": r.mean(), "COD": cod, "PRD": prd, "PRB": b[1], "PRB_se": se,77            "ratio_pondere": V.sum() / P.sum()}7879# ----------------------------------------------------------------------------- spatial8081def matrice_poids_knn(coords: np.ndarray, k: int = 8, standardiser: bool = True) -> np.ndarray:82    """Matrice W (n×n) des k plus proches voisins (distance euclidienne), standardisée en ligne."""83    c = np.asarray(coords, float)84    d = np.sqrt(((c[:, None, :] - c[None, :, :]) ** 2).sum(-1))85    np.fill_diagonal(d, np.inf)86    W = np.zeros_like(d)87    idx = np.argsort(d, axis=1)[:, :k]88    for i in range(len(c)):89        W[i, idx[i]] = 1.090    if standardiser:91        W = W / W.sum(1, keepdims=True)92    return W9394def matrice_poids_distance_inverse(coords: np.ndarray, portee: float, puissance: float = 1.0, standardiser: bool = True) -> np.ndarray:95    c = np.asarray(coords, float)96    d = np.sqrt(((c[:, None, :] - c[None, :, :]) ** 2).sum(-1))97    with np.errstate(divide="ignore"):98        W = np.where((d > 0) & (d <= portee), 1.0 / d ** puissance, 0.0)99    if standardiser:100        s = W.sum(1, keepdims=True); s[s == 0] = 1101        W = W / s102    return W103104def moran_i(x: np.ndarray, W: np.ndarray) -> dict:105    """I de Moran et test sous hypothèse de normalité : I = (n/S0) · (z'Wz)/(z'z)."""106    z = np.asarray(x, float); z = z - z.mean()107    n = len(z); S0 = W.sum()108    I = (n / S0) * (z @ W @ z) / (z @ z)109    EI = -1.0 / (n - 1)110    S1 = 0.5 * ((W + W.T) ** 2).sum()111    S2 = ((W.sum(0) + W.sum(1)) ** 2).sum()112    VI = (n ** 2 * S1 - n * S2 + 3 * S0 ** 2) / ((n ** 2 - 1) * S0 ** 2) - EI ** 2113    zI = (I - EI) / np.sqrt(VI)114    from math import erf, sqrt115    p = 2 * (1 - 0.5 * (1 + erf(abs(zI) / sqrt(2))))116    return {"I": I, "E_I": EI, "Var_I": VI, "z": zI, "p": p}117118def decalage_spatial(W: np.ndarray, x: np.ndarray) -> np.ndarray:119    return W @ np.asarray(x, float)120121def sar_lag_mv(y: np.ndarray, X: np.ndarray, W: np.ndarray, grille_rho=None) -> dict:122    """Modèle autorégressif spatial y = ρWy + Xβ + ε estimé par maximum de vraisemblance concentrée (recherche sur ρ)."""123    y = np.asarray(y, float); X = np.asarray(X, float); n = len(y)124    Wy = W @ y125    b0, *_ = np.linalg.lstsq(X, y, rcond=None); e0 = y - X @ b0126    bL, *_ = np.linalg.lstsq(X, Wy, rcond=None); eL = Wy - X @ bL127    eig = np.linalg.eigvals(W).real128    rhos = grille_rho if grille_rho is not None else np.linspace(-0.9, 0.99, 380)129    best = None130    for rho in rhos:131        e = e0 - rho * eL132        s2 = e @ e / n133        ll = np.log(np.abs(1 - rho * eig)).sum() - n / 2 * np.log(2 * np.pi * s2) - n / 2134        if best is None or ll > best[0]:135            best = (ll, rho, s2)136    ll, rho, s2 = best137    beta = b0 - rho * bL138    return {"rho": rho, "beta": beta, "sigma2": s2, "loglik": ll}139