#!/usr/bin/env python3
"""Chapitre 18 — Modèles hédoniques : exemple intégral sur les 600 ventes du cas 1 (Gatineau, DONNÉES SYNTHÉTIQUES, graine 20261031).

Le modèle de référence est `traite.cas1.hedonique` (semi-log, HC1) ; il n'est PAS redéfini ici (les ch. 17, 19, 20, 35, 37, 39, 45 s'y réfèrent).
Ce script ajoute : la comparaison des coefficients estimés aux paramètres de génération (data/synthetic/README.md), les variantes de
spécification (linéaire, log-log, Box–Cox, sans secteurs, indicatrices trimestrielles), les diagnostics (Breusch–Pagan, Jarque–Bera, VIF,
distances de Cook, I de Moran des résidus), les prix implicites au sujet (en $), la prédiction avec correction de Duan et son intervalle,
les statistiques de ratio (COD, PRD, PRB) du modèle et du rôle fictif, et l'analyse de sensibilité (spécification, sous-échantillons).
Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch18_hedonique.py"""
from __future__ import annotations
import sys, math, pathlib
import numpy as np, pandas as pd
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
from traite import cas1
from traite import hedonique as h
from scipy import stats, optimize

ROOT = pathlib.Path(__file__).resolve().parents[3]
GRAINE = 20261031

# ----------------------------------------------------------------------------- paramètres de génération (README synthétique)
_g = (1 + 9.06 / 100) ** 0.5 - 1          # croissance annuelle composée du jeu (variation du rôle 2027-2029 sur deux ans)
VRAIS = {"ln_S": 0.55, "ln_T": 0.10, "age": -0.0055, "age2_100": 0.0030, "qualite": 0.06, "garage": 0.035, "sous_sol": 0.03,
         "piscine": 0.02, "etages": 0.02, "dist_centre": -0.006,
         "sect_Gatineau": math.log(464_900 / 475_000), "sect_Aylmer": math.log(530_000 / 475_000),
         "sect_Buckingham": math.log(395_000 / 475_000), "sect_Masson-Angers": math.log(395_000 / 475_000) - 0.02,
         "mois": math.log1p(_g) / 12, "sigma": 0.085}
ETIQUETTES = {"const": "Constante", "ln_S": "ln(superficie habitable / 140 m²)", "ln_T": "ln(terrain / 560 m²)", "age": "Âge (années)",
              "age2_100": "Âge² / 100", "qualite": "Qualité (écart à 3)", "garage": "Garage (places)", "sous_sol": "Sous-sol fini",
              "piscine": "Piscine", "etages": "À étages", "dist_centre": "Distance au centre au-delà de 3 km (km)",
              "sect_Gatineau": "Secteur Gatineau", "sect_Aylmer": "Secteur Aylmer", "sect_Buckingham": "Secteur Buckingham",
              "sect_Masson-Angers": "Secteur Masson-Angers", "mois": "Mois depuis le 1er juillet 2024"}

# ----------------------------------------------------------------------------- utilitaires
def _duan(res: np.ndarray) -> float:
    return float(np.mean(np.exp(res)))

def ratios(valeurs: np.ndarray, prix: np.ndarray) -> dict:
    return h.statistiques_ratio(valeurs, prix)

def _rmse(a, b) -> float:
    a = np.asarray(a, float); b = np.asarray(b, float)
    return float(np.sqrt(np.mean((a - b) ** 2)))

def trimestre(d: pd.Series) -> np.ndarray:
    d = pd.to_datetime(d)
    return ((d.dt.year - 2023) * 4 + (d.dt.month - 1) // 3).to_numpy()

# ----------------------------------------------------------------------------- variantes de spécification
def matrice_lineaire(v: pd.DataFrame) -> pd.DataFrame:
    """Forme linéaire : mêmes variables, superficie et terrain en niveaux (m²)."""
    X = cas1.matrice(v).copy()
    X["S"] = v["superficie_habitable_m2"].astype(float); X["T"] = v["superficie_terrain_m2"].astype(float)
    X = X.drop(columns=["ln_S", "ln_T"])
    return X[["const", "S", "T"] + [c for c in X.columns if c not in ("const", "S", "T")]]

def matrice_loglog(v: pd.DataFrame) -> pd.DataFrame:
    """Forme log-log : l'âge entre en logarithme (ln(1 + âge)) à la place du polynôme."""
    X = cas1.matrice(v).copy()
    X["ln_age"] = np.log1p(v["age"].astype(float))
    X = X.drop(columns=["age", "age2_100"])
    return X[["const", "ln_S", "ln_T", "ln_age"] + [c for c in X.columns if c not in ("const", "ln_S", "ln_T", "ln_age")]]

def matrice_sans_secteurs(v: pd.DataFrame) -> pd.DataFrame:
    X = cas1.matrice(v)
    return X[[c for c in X.columns if not c.startswith("sect_")]]

def matrice_trimestres(v: pd.DataFrame) -> pd.DataFrame:
    """Indicatrices trimestrielles (2023T1 = référence) à la place de la tendance linéaire `mois`."""
    X = cas1.matrice(v).drop(columns=["mois"])
    q = trimestre(v["date_vente"])
    for k in range(1, int(q.max()) + 1):
        X[f"trim_{k}"] = (q == k).astype(float)
    return X

def boxcox_profil(P: np.ndarray, X: np.ndarray) -> dict:
    """Transformation de Box–Cox (1964) de la variable dépendante : y(λ) = (P^λ − 1)/λ, λ = 0 → ln P.
    Log-vraisemblance concentrée : ℓ(λ) = −(n/2) ln(SCR_λ / n) + (λ − 1) Σ ln P. IC à 95 % : ℓ(λ) ≥ ℓ_max − χ²₁(0,95)/2."""
    P = np.asarray(P, float); n = len(P); slog = np.log(P).sum()
    def y_l(l): return np.log(P) if abs(l) < 1e-8 else (P ** l - 1) / l
    def ll(l):
        y = y_l(l); b = np.linalg.lstsq(X, y, rcond=None)[0]; r = y - X @ b
        return -(n / 2) * np.log(r @ r / n) + (l - 1) * slog
    opt = optimize.minimize_scalar(lambda l: -ll(l), bounds=(-1.0, 1.5), method="bounded")
    lam = float(opt.x); llmax = ll(lam); seuil = llmax - stats.chi2.ppf(0.95, 1) / 2
    bas = optimize.brentq(lambda l: ll(l) - seuil, -1.0, lam); haut = optimize.brentq(lambda l: ll(l) - seuil, lam, 1.5)
    grille = np.linspace(-0.6, 1.2, 73)
    return {"lambda": lam, "ic95": (float(bas), float(haut)), "ll_max": float(llmax), "ll_0": float(ll(0.0)), "ll_1": float(ll(1.0)),
            "grille": grille, "ll_grille": np.array([ll(l) for l in grille]), "y": y_l(lam)}

def _predire_lineaire(beta, X): return X @ beta

def ajuster_variantes(v: pd.DataFrame, s: dict, ref: dict) -> dict:
    """Estime les variantes et renvoie, pour chacune : R², σ, RMSE en $, COD, prédiction du sujet."""
    P = v["prix_vente"].astype(float).to_numpy(); lnP = np.log(P)
    sd = cas1.sujet_df(s)
    out = {}
    # --- référence (semi-log)
    out["semi-log (référence)"] = {"r2": ref["r2"], "r2_ajuste": ref["r2_ajuste"], "sigma": ref["sigma"], "k": ref["k"],
                                   "ajuste": np.exp(ref["X"].values @ ref["beta"]) * _duan(ref["residus"]), "sujet": float(cas1.predire(ref, sd)[0])}
    # --- linéaire
    XL = matrice_lineaire(v); rL = h.mco(XL.values, P)
    XLs = matrice_lineaire(sd)[XL.columns].values
    out["linéaire"] = {"r2": rL["r2"], "r2_ajuste": rL["r2_ajuste"], "sigma": float(np.sqrt(rL["sigma2"])), "k": rL["k"],
                       "ajuste": XL.values @ rL["beta"], "sujet": float((XLs @ rL["beta"])[0]), "mco": rL, "X": XL}
    # --- log-log
    XG = matrice_loglog(v); rG = h.mco(XG.values, lnP)
    XGs = matrice_loglog(sd)[XG.columns].values
    out["log-log"] = {"r2": rG["r2"], "r2_ajuste": rG["r2_ajuste"], "sigma": float(np.sqrt(rG["sigma2"])), "k": rG["k"],
                      "ajuste": np.exp(XG.values @ rG["beta"]) * _duan(rG["residus"]), "sujet": float(np.exp(XGs @ rG["beta"])[0] * _duan(rG["residus"])),
                      "mco": rG, "X": XG}
    # --- Box–Cox
    X0 = ref["X"].values; bc = boxcox_profil(P, X0)
    rB = h.mco(X0, bc["y"]); lam = bc["lambda"]
    def inv(yl): return np.exp(yl) if abs(lam) < 1e-8 else np.power(np.maximum(lam * yl + 1, 1e-9), 1 / lam)
    Xs0 = cas1.matrice(sd)[ref["X"].columns].values
    # retransformation par lissage de Duan généralisé : moyenne des inv(ŷ + e_i)
    aj_bc = np.array([np.mean(inv(yh + rB["residus"])) for yh in (X0 @ rB["beta"])])
    out["Box–Cox"] = {"r2": rB["r2"], "r2_ajuste": rB["r2_ajuste"], "sigma": float(np.sqrt(rB["sigma2"])), "k": rB["k"], "lambda": lam,
                      "ic95": bc["ic95"], "ajuste": aj_bc, "sujet": float(np.mean(inv((Xs0 @ rB["beta"])[0] + rB["residus"]))), "mco": rB, "boxcox": bc}
    # --- sans secteurs
    XS = matrice_sans_secteurs(v); rS = h.mco(XS.values, lnP)
    XSs = matrice_sans_secteurs(sd)[XS.columns].values
    out["semi-log sans secteurs"] = {"r2": rS["r2"], "r2_ajuste": rS["r2_ajuste"], "sigma": float(np.sqrt(rS["sigma2"])), "k": rS["k"],
                                     "ajuste": np.exp(XS.values @ rS["beta"]) * _duan(rS["residus"]), "sujet": float(np.exp(XSs @ rS["beta"])[0] * _duan(rS["residus"])),
                                     "mco": rS, "X": XS, "coef": pd.Series(rS["beta"], index=XS.columns), "se": pd.Series(rS["se_hc1"], index=XS.columns)}
    # --- indicatrices trimestrielles
    XT = matrice_trimestres(v); rT = h.mco(XT.values, lnP)
    # sujet : date d'évaluation 2026-09-01 → trimestre 2026T3 (= dernier trimestre observé)
    sdT = cas1.matrice(sd).drop(columns=["mois"])
    for c in XT.columns:
        if c.startswith("trim_"): sdT[c] = 0.0
    sdT[f"trim_{int(trimestre(sd['date_vente']).max())}"] = 1.0
    XTs = sdT[XT.columns].values
    out["semi-log, indicatrices trimestrielles"] = {"r2": rT["r2"], "r2_ajuste": rT["r2_ajuste"], "sigma": float(np.sqrt(rT["sigma2"])), "k": rT["k"],
                                                   "ajuste": np.exp(XT.values @ rT["beta"]) * _duan(rT["residus"]), "sujet": float(np.exp(XTs @ rT["beta"])[0] * _duan(rT["residus"])),
                                                   "mco": rT, "X": XT}
    for nom, d in out.items():
        d["rmse"] = _rmse(d["ajuste"], P); d["ratios"] = ratios(d["ajuste"], P)
    return out

# ----------------------------------------------------------------------------- diagnostics
def breusch_pagan(X: np.ndarray, res: np.ndarray) -> dict:
    """Test de Breusch–Pagan (1979), version studentisée de Koenker : LM = n·R² de la régression de e² sur X ; χ²(k − 1)."""
    e2 = res ** 2
    r = h.mco(X, e2, robuste=False)
    n, k = X.shape; LM = n * r["r2"]
    return {"LM": float(LM), "ddl": k - 1, "p": float(1 - stats.chi2.cdf(LM, k - 1))}

def vif(X: pd.DataFrame) -> pd.Series:
    cols = [c for c in X.columns if c != "const"]
    out = {}
    for c in cols:
        autres = X[[k for k in X.columns if k != c]].values
        r = h.mco(autres, X[c].values, robuste=False)
        out[c] = 1.0 / max(1.0 - r["r2"], 1e-12)
    return pd.Series(out)

def cook(X: np.ndarray, res: np.ndarray, sigma2: float) -> dict:
    n, k = X.shape
    H = X @ np.linalg.pinv(X.T @ X) @ X.T
    hii = np.diag(H)
    D = res ** 2 * hii / (k * sigma2 * (1 - hii) ** 2)
    return {"D": D, "h": hii, "seuil": 4.0 / n, "n_sup": int((D > 4.0 / n).sum())}

def moran_residus(v: pd.DataFrame, res: np.ndarray, k: int = 8) -> dict:
    coords = np.column_stack([(v["lat"] - 45.428) * 111.0, (v["lon"] + 75.714) * 78.0])
    W = h.matrice_poids_knn(coords, k)
    return h.moran_i(res, W)

def diagnostics(v: pd.DataFrame, ref: dict) -> dict:
    X = ref["X"]; res = ref["residus"]
    jb = stats.jarque_bera(res)
    return {"breusch_pagan": breusch_pagan(X.values, res), "jarque_bera": {"stat": float(jb.statistic), "p": float(jb.pvalue)},
            "skew": float(stats.skew(res)), "kurtosis": float(stats.kurtosis(res, fisher=False)),
            "vif": vif(X), "cook": cook(X.values, res, ref["sigma"] ** 2), "moran": moran_residus(v, res),
            "res_min": float(res.min()), "res_max": float(res.max()), "pct_abs_sup_15": float(np.mean(np.abs(res) > 0.15) * 100)}

# ----------------------------------------------------------------------------- prix implicites au sujet
def prix_implicites_sujet(ref: dict, s: dict) -> pd.DataFrame:
    """Prix implicites en dollars au sujet, dérivés du semi-log autour de la prédiction corrigée P̂_s (Duan).
    Continues : ∂P/∂x = β P̂ (× dérivée de la transformation) ; indicatrices : P̂ [exp(β) − 1] (Halvorsen–Palmquist) et Kennedy."""
    beta = pd.Series(ref["beta"], index=ref["X"].columns); se = pd.Series(np.sqrt(np.diag(ref["V_hc1"])), index=ref["X"].columns)
    sd = cas1.sujet_df(s); Ps = float(cas1.predire(ref, sd)[0]); xs = cas1.matrice(sd).iloc[0]
    S = float(s["superficie_habitable_m2"]); T = float(s["superficie_terrain_m2"]); age = float(s["age"])
    rows = []
    def cont(nom, unite, b, scale, se_b):
        rows.append({"caracteristique": nom, "unite": unite, "type": "continue", "beta": b, "prix_implicite": b * Ps * scale, "se": abs(scale) * Ps * se_b})
    cont("Superficie habitable", "$/m²", beta["ln_S"], 1.0 / S, se["ln_S"])
    cont("Superficie du terrain", "$/m²", beta["ln_T"], 1.0 / T, se["ln_T"])
    d_age = beta["age"] + 2 * beta["age2_100"] * age / 100.0
    se_age = math.sqrt(ref["V_hc1"][list(beta.index).index("age"), list(beta.index).index("age")]
                       + (2 * age / 100) ** 2 * ref["V_hc1"][list(beta.index).index("age2_100"), list(beta.index).index("age2_100")]
                       + 2 * (2 * age / 100) * ref["V_hc1"][list(beta.index).index("age"), list(beta.index).index("age2_100")])
    rows.append({"caracteristique": "Âge (une année de plus, à 28 ans)", "unite": "$/an", "type": "continue", "beta": d_age, "prix_implicite": d_age * Ps, "se": Ps * se_age})
    cont("Conditions du marché (un mois)", "$/mois", beta["mois"], 1.0, se["mois"])
    rows.append({"caracteristique": "Distance au centre (sujet à 2,1 km : dans le plateau)", "unite": "$/km", "type": "continue", "beta": beta["dist_centre"], "prix_implicite": 0.0, "se": 0.0})
    def disc(nom, col, direction=+1):
        b = beta[col] * direction; v = se[col] ** 2
        hp = Ps * (math.exp(b) - 1); ken = Ps * (math.exp(b - v / 2) - 1)
        rows.append({"caracteristique": nom, "unite": "$", "type": "indicatrice", "beta": beta[col], "prix_implicite": hp, "kennedy": ken, "se": Ps * math.exp(b) * se[col]})
    rows.append({"caracteristique": "Qualité (un échelon de plus)", "unite": "$", "type": "indicatrice", "beta": beta["qualite"],
                 "prix_implicite": Ps * (math.exp(beta["qualite"]) - 1), "kennedy": Ps * (math.exp(beta["qualite"] - se["qualite"] ** 2 / 2) - 1), "se": Ps * math.exp(beta["qualite"]) * se["qualite"]})
    disc("Garage (une place de plus)", "garage")
    disc("Sous-sol fini (contribution : fini vs non fini)", "sous_sol")
    disc("Piscine (ajout)", "piscine")
    disc("Cottage à étages vs plain-pied", "etages")
    disc("Secteur Aylmer vs Hull", "sect_Aylmer")
    disc("Secteur Gatineau vs Hull", "sect_Gatineau")
    df = pd.DataFrame(rows); df["P_sujet"] = Ps
    return df

# ----------------------------------------------------------------------------- prédiction et intervalle
def prediction_sujet(ref: dict, s: dict) -> dict:
    sd = cas1.sujet_df(s); X = cas1.matrice(sd)[ref["X"].columns].values
    lnp = float((X @ ref["beta"])[0]); xVx = float((X @ ref["V_hc1"] @ X.T)[0, 0])
    duan = _duan(ref["residus"]); logn = math.exp(ref["sigma"] ** 2 / 2)
    ip = cas1.intervalle_prediction(ref, sd).iloc[0]
    ip80 = cas1.intervalle_prediction(ref, sd, 0.80).iloc[0]
    # prédiction à la date de référence du rôle fictif (2024-07-01) pour la comparaison avec la valeur au rôle
    sd_ref = cas1.sujet_df(s, "2024-07-01"); p_ref = float(cas1.predire(ref, sd_ref)[0])
    return {"ln_pred": lnp, "naif": math.exp(lnp), "duan": math.exp(lnp) * duan, "facteur_duan": duan, "facteur_lognormal": logn,
            "lognormal": math.exp(lnp) * logn, "xVx": xVx, "sigma2": ref["sigma"] ** 2, "sd_pred": math.sqrt(ref["sigma"] ** 2 + xVx),
            "ic95": (float(ip["bas"]), float(ip["haut"])), "ic80": (float(ip80["bas"]), float(ip80["haut"])),
            "pred_2024_07": p_ref, "ratio_role": float(s["valeur_role_2025"]) / p_ref}

# ----------------------------------------------------------------------------- validation croisée et ratios
def validation_croisee(v: pd.DataFrame, nfolds: int = 10, graine: int = GRAINE) -> dict:
    rng = np.random.default_rng(graine)
    X = cas1.matrice(v).values; lnP = np.log(v["prix_vente"].astype(float).to_numpy()); P = np.exp(lnP)
    idx = rng.permutation(len(P)); folds = np.array_split(idx, nfolds)
    pred = np.empty_like(P)
    for f in folds:
        tr = np.setdiff1d(idx, f); r = h.mco(X[tr], lnP[tr])
        pred[f] = np.exp(X[f] @ r["beta"]) * _duan(r["residus"])
    return {"pred": pred, "rmse": _rmse(pred, P), "ratios": ratios(pred, P), "erreur_abs_pct_mediane": float(np.median(np.abs(pred / P - 1)) * 100)}

def ratios_role_fictif(v: pd.DataFrame, ref: dict) -> dict:
    """Étude de ratios du rôle fictif : valeur au rôle / prix ramené à la date de référence (2024-07-01) par la tendance du modèle."""
    bm = ref["beta"][list(ref["X"].columns).index("mois")]
    prix_ref = v["prix_vente"].astype(float).to_numpy() * np.exp(-bm * ref["X"]["mois"].to_numpy())
    return {"ratios": ratios(v["valeur_role_2025"].to_numpy(float), prix_ref), "prix_ref": prix_ref,
            "r": v["valeur_role_2025"].to_numpy(float) / prix_ref}

def indice_trimestriel(variantes: dict, ref: dict, v: pd.DataFrame) -> pd.DataFrame:
    rT = variantes["semi-log, indicatrices trimestrielles"]["mco"]; XT = variantes["semi-log, indicatrices trimestrielles"]["X"]
    cols = [c for c in XT.columns if c.startswith("trim_")]
    bq = np.array([rT["beta"][list(XT.columns).index(c)] for c in cols]); seq = np.array([rT["se_hc1"][list(XT.columns).index(c)] for c in cols])
    idx = h.indice_dummies_temps(bq)
    q = np.arange(len(cols) + 1)
    # tendance : mois au milieu de chaque trimestre depuis 2024-07-01, normalisée à 2023T1
    mid = pd.to_datetime([f"{2023 + k // 4}-{3 * (k % 4) + 2:02d}-15" for k in q])
    mois_mid = (mid - cas1.REF_DATE).days / 30.4375
    bm = ref["beta"][list(ref["X"].columns).index("mois")]
    tend = 100 * np.exp(bm * (mois_mid - mois_mid[0]))
    lab = [f"{2023 + k // 4}T{k % 4 + 1}" for k in q]
    return pd.DataFrame({"trimestre": lab, "indice_dummies": idx, "bas": np.concatenate([[100], 100 * np.exp(bq - 1.96 * seq)]),
                         "haut": np.concatenate([[100], 100 * np.exp(bq + 1.96 * seq)]), "indice_tendance": tend,
                         "n_ventes": [int((trimestre(v["date_vente"]) == k).sum()) for k in q]})

# ----------------------------------------------------------------------------- sensibilité : sous-échantillons
def sensibilite_echantillon(v: pd.DataFrame, s: dict, ref: dict, diag: dict) -> pd.DataFrame:
    sd = cas1.sujet_df(s); Xs_full = cas1.matrice(sd)
    def fit(sub: pd.DataFrame, nom: str):
        X = cas1.matrice(sub); y = np.log(sub["prix_vente"].astype(float)).values
        garder = [c for c in X.columns if X[c].std() > 0 or c == "const"]
        X = X[garder]; r = h.mco(X.values, y)
        p = float(np.exp((Xs_full[garder].values @ r["beta"])[0]) * _duan(r["residus"]))
        bS = r["beta"][garder.index("ln_S")]
        return {"echantillon": nom, "n": len(sub), "r2": r["r2"], "sigma": float(np.sqrt(r["sigma2"])), "beta_ln_S": bS, "sujet": p}
    rows = [fit(v, "600 ventes (référence)")]
    rows.append(fit(v[v["secteur"] == "Hull"], "Secteur Hull seulement"))
    mois_av = (pd.Timestamp(s["date_evaluation"]) - v["date_vente"]).dt.days / 30.4375
    rows.append(fit(v[mois_av <= 18], "Ventes des 18 derniers mois"))
    rows.append(fit(v[v["type"] == "à étages"], "Cottages à étages seulement"))
    rows.append(fit(v[diag["cook"]["D"] <= diag["cook"]["seuil"]], "Sans les observations influentes (Cook > 4/n)"))
    rows.append(fit(v[(v["superficie_habitable_m2"].between(110, 230))], "Superficie de 110 à 230 m²"))
    return pd.DataFrame(rows)

# ----------------------------------------------------------------------------- exercices (chiffres des solutions)
def exercices(ref: dict, pred: dict, vari: dict, idx: pd.DataFrame) -> dict:
    ex = {}
    # 18.1 : interprétation des coefficients
    ex["1"] = {"elasticite_S": 0.539, "garage_hp": math.exp(0.0324) - 1, "buck_hp": math.exp(-0.250) - 1, "buck_kennedy": math.exp(-0.250 - 0.072 ** 2 / 2) - 1,
               "mois_annuel": (1 + 0.0036) ** 12 - 1}
    # 18.2 : prix implicites à partir de coefficients arrondis
    P, S = 638_000.0, 165.0
    ex["2"] = {"pS_marginal": 0.539 * P / S, "pS_10m2_exact": P * ((175 / 165) ** 0.539 - 1), "pS_10m2_lineaire": 0.539 * P / S * 10,
               "age28": (-0.0051 + 2 * 0.0026 * 28 / 100) * P, "age60": (-0.0051 + 2 * 0.0026 * 60 / 100) * P, "age_min": 0.0051 / (2 * 0.0026 / 100),
               "piscine_hp": P * (math.exp(0.0263) - 1)}
    # 18.3 : Box–Cox
    ex["3"] = {"bc_1": (500_000 ** 1 - 1) / 1, "bc_05": (500_000 ** 0.5 - 1) / 0.5, "bc_0": math.log(500_000), "lambda_hat": vari["Box–Cox"]["lambda"], "ic": vari["Box–Cox"]["ic95"]}
    # 18.4 : Duan et intervalle avec chiffres arrondis
    lnp, sig, xVx, duan = 13.3627, 0.0830, round(pred["xVx"], 5), 1.0034
    sd = math.sqrt(sig ** 2 + xVx); z = 1.959964
    ex["4"] = {"naif": math.exp(lnp), "duan": math.exp(lnp) * duan, "lognormal": math.exp(lnp) * math.exp(sig ** 2 / 2), "sd": sd,
               "bas": math.exp(lnp - z * sd), "haut": math.exp(lnp + z * sd), "xVx": xVx, "part_xVx": xVx / (sig ** 2 + xVx),
               "demi_largeur_pct": (math.exp(z * sd) - 1) * 100}
    # 18.5 : étude de ratios (dix ventes fictives)
    prix = np.array([312_000, 358_500, 401_000, 445_000, 489_500, 532_000, 598_000, 655_000, 742_000, 895_000], float)
    val = np.array([298_000, 329_000, 371_000, 405_000, 452_000, 478_000, 539_000, 576_000, 645_000, 752_000], float)
    r5 = ratios(val, prix); ex["5"] = {"ratios": val / prix, **r5}
    # 18.6 : variables omises (coefficients sans secteurs)
    cS = vari["semi-log sans secteurs"]["coef"]; cR = pd.Series(ref["beta"], index=ref["X"].columns)
    ex["6"] = {"dist_ref": cR["dist_centre"], "dist_sans": cS["dist_centre"], "lnT_ref": cR["ln_T"], "lnT_sans": cS["ln_T"],
               "sigma_ref": ref["sigma"], "sigma_sans": vari["semi-log sans secteurs"]["sigma"], "sujet_ref": vari["semi-log (référence)"]["sujet"],
               "sujet_sans": vari["semi-log sans secteurs"]["sujet"], "r2_sans": vari["semi-log sans secteurs"]["r2"]}
    # variante demandée à l'exercice : omission de la qualité (tirée indépendamment des autres attributs dans la loi de génération)
    XQ = ref["X"].drop(columns=["qualite"]); rQ = h.mco(XQ.values, ref["y"])
    d = cas1.charger(); sdQ = cas1.matrice(cas1.sujet_df(d["sujet"]))[XQ.columns].values
    ex["6"].update({"sigma_sans_qualite": float(np.sqrt(rQ["sigma2"])), "r2_sans_qualite": rQ["r2"], "lnS_sans_qualite": rQ["beta"][list(XQ.columns).index("ln_S")],
                    "sujet_sans_qualite": float(np.exp(sdQ @ rQ["beta"])[0] * _duan(rQ["residus"]))})
    # 18.7 : Rosen — prix implicite marginal et enchère
    ex["7"] = {"p_S": 0.539 * 638_000 / 165, "p_S_120": 0.539 * 638_000 * (120 / 165) ** 0.539 / 120, "p_S_220": 0.539 * 638_000 * (220 / 165) ** 0.539 / 220}
    # 18.8 : indice trimestriel vs tendance
    bm = ref["beta"][list(ref["X"].columns).index("mois")]
    ecart = (idx["indice_dummies"] / idx["indice_tendance"] - 1) * 100
    ex["8"] = {"g_mensuel": math.exp(bm) - 1, "g_annuel": math.exp(12 * bm) - 1, "ecart_abs_moyen_pct": float(ecart.abs().mean()), "ecart_max_pct": float(ecart.abs().max()),
               "G0404_tendance": 506_500 * (math.exp(bm * 15.2) - 1), "indice_fin": float(idx["indice_dummies"].iloc[-1]), "tendance_fin": float(idx["indice_tendance"].iloc[-1]),
               "indice_2025T2": float(idx.set_index("trimestre").loc["2025T2", "indice_dummies"]),
               "G0404_indices": 506_500 * (float(idx["indice_dummies"].iloc[-1]) / float(idx.set_index("trimestre").loc["2025T2", "indice_dummies"]) - 1)}
    return ex

# ----------------------------------------------------------------------------- principal
def main(verbose: bool = True) -> dict:
    d = cas1.charger(); v = d["ventes"]; s = d["sujet"]
    ref = cas1.hedonique(v)
    coef = ref["coef"].copy(); coef["vrai"] = [VRAIS.get(c, np.nan) for c in coef.index]
    coef["ecart_t"] = (coef["coef"] - coef["vrai"]) / coef["se_hc1"]
    coef["ic_bas"] = coef["coef"] - 1.96 * coef["se_hc1"]; coef["ic_haut"] = coef["coef"] + 1.96 * coef["se_hc1"]
    coef["couvre"] = (coef["ic_bas"] <= coef["vrai"]) & (coef["vrai"] <= coef["ic_haut"])
    coef["etiquette"] = [ETIQUETTES[c] for c in coef.index]
    diag = diagnostics(v, ref)
    vari = ajuster_variantes(v, s, ref)
    pi = prix_implicites_sujet(ref, s)
    pred = prediction_sujet(ref, s)
    cv = validation_croisee(v)
    role = ratios_role_fictif(v, ref)
    idx = indice_trimestriel(vari, ref, v)
    sens = sensibilite_echantillon(v, s, ref, diag)
    ex = exercices(ref, pred, vari, idx)
    # comparables du ch. 16 : résidus du modèle
    sel = cas1.selection(); comps = sel["comparables"]
    res_comps = pd.Series(ref["residus"], index=v.index).loc[comps.index]
    P = v["prix_vente"].astype(float).to_numpy()
    out = {"ventes": v, "sujet": s, "modele": ref, "coef": coef, "diag": diag, "variantes": vari, "prix_implicites": pi, "prediction": pred,
           "cv": cv, "role": role, "indice": idx, "sensibilite": sens, "exercices": ex, "comparables": comps, "residus_comparables": res_comps,
           "ratios_modele": ratios(vari["semi-log (référence)"]["ajuste"], P)}
    if verbose:
        pd.set_option("display.width", 220); pd.set_option("display.max_columns", 30)
        print("== Modèle de référence (traite.cas1.hedonique) : n =", ref["n"], "k =", ref["k"], f"R² = {ref['r2']:.4f}  R²aj = {ref['r2_ajuste']:.4f}  σ = {ref['sigma']:.4f}")
        print(coef[["coef", "se_hc1", "t", "effet_pct", "vrai", "ecart_t", "couvre"]].round(4).to_string())
        print(f"   couverture des IC 95 % : {int(coef['couvre'].sum())}/{int(coef['vrai'].notna().sum())} paramètres de génération couverts")
        print(f"== Croissance mensuelle {ref['croissance_mensuelle']*100:.3f} % ; annuelle {((1+ref['croissance_mensuelle'])**12-1)*100:.2f} % ; vraie {_g*100:.2f} %")
        bp = diag["breusch_pagan"]; print(f"== Breusch–Pagan (Koenker) LM = {bp['LM']:.2f}, ddl = {bp['ddl']}, p = {bp['p']:.3f}")
        print(f"== Jarque–Bera = {diag['jarque_bera']['stat']:.2f} (p = {diag['jarque_bera']['p']:.3f}) ; asymétrie {diag['skew']:.3f} ; aplatissement {diag['kurtosis']:.3f}")
        print(f"== Résidus : min {diag['res_min']:.3f}, max {diag['res_max']:.3f}, |e| > 0,15 : {diag['pct_abs_sup_15']:.1f} %")
        print("== VIF :\n", diag["vif"].round(2).to_string())
        ck = diag["cook"]; top = np.argsort(-ck["D"])[:3]
        print(f"== Cook : seuil 4/n = {ck['seuil']:.4f} ; {ck['n_sup']} observations au-dessus ; max {ck['D'].max():.4f} ; 3 plus grandes :", list(v['id_vente'].iloc[top]), np.round(ck['D'][top], 4))
        m = diag["moran"]; print(f"== Moran (résidus, 8 voisins) I = {m['I']:.4f}, z = {m['z']:.2f}, p = {m['p']:.3f}")
        print("== Variantes :")
        for nom, dd in vari.items():
            extra = f" λ = {dd['lambda']:.3f} IC [{dd['ic95'][0]:.3f}; {dd['ic95'][1]:.3f}]" if "lambda" in dd else ""
            print(f"   {nom:42s} k={dd['k']:2d} R²={dd['r2']:.4f} R²aj={dd['r2_ajuste']:.4f} σ={dd['sigma']:.4f} RMSE={dd['rmse']:,.0f} $ COD={dd['ratios']['COD']:.2f} PRD={dd['ratios']['PRD']:.3f} PRB={dd['ratios']['PRB']:.3f} sujet={dd['sujet']:,.0f} ${extra}")
        cS = vari["semi-log sans secteurs"]
        print(f"   sans secteurs : dist_centre {cS['coef']['dist_centre']:.4f} (réf {ref['beta'][list(ref['X'].columns).index('dist_centre')]:.4f}) ; ln_T {cS['coef']['ln_T']:.4f} (réf {ref['beta'][list(ref['X'].columns).index('ln_T')]:.4f})")
        print("== Prix implicites au sujet (P̂ = %.0f $) :" % pi["P_sujet"].iloc[0])
        print(pi[["caracteristique", "unite", "beta", "prix_implicite", "kennedy", "se"]].round(4).to_string(index=False))
        print("== Prédiction du sujet :", {k: (round(x, 5) if isinstance(x, float) and abs(x) < 10 else (tuple(round(y) for y in x) if isinstance(x, tuple) else round(x))) for k, x in pred.items()})
        print(f"== Ratios du modèle (in-sample) : COD {out['ratios_modele']['COD']:.2f}, PRD {out['ratios_modele']['PRD']:.4f}, PRB {out['ratios_modele']['PRB']:.4f}, médiane {out['ratios_modele']['ratio_median']:.4f}")
        print(f"== Validation croisée 10 plis : RMSE {cv['rmse']:,.0f} $, COD {cv['ratios']['COD']:.2f}, PRD {cv['ratios']['PRD']:.4f}, PRB {cv['ratios']['PRB']:.4f}, erreur absolue médiane {cv['erreur_abs_pct_mediane']:.2f} %")
        rr = role["ratios"]; print(f"== Rôle fictif vs prix ramenés au 2024-07-01 : médiane {rr['ratio_median']:.4f}, COD {rr['COD']:.2f}, PRD {rr['PRD']:.4f}, PRB {rr['PRB']:.4f}")
        print("== Indice trimestriel :\n", idx.round(2).to_string(index=False))
        print("== Sensibilité à l'échantillon :\n", sens.round(4).to_string(index=False))
        print("== Résidus des 5 comparables du ch. 16 :", dict(zip(comps["id_vente"], res_comps.round(4))), "moyenne", round(res_comps.mean(), 4))
        print("== Exercices :")
        for k, e in ex.items(): print("  ", k, {kk: (np.round(vv, 4) if not isinstance(vv, tuple) else tuple(round(x, 4) for x in vv)) for kk, vv in e.items()})
    return out

if __name__ == "__main__":
    main()