#!/usr/bin/env python3
"""Chapitre 35 — Indices immobiliers et ventes répétées : exemple intégral sur le cas 1 (Gatineau, DONNÉES SYNTHÉTIQUES, graine 20261031).

Construit et compare, sur les 600 ventes du cas 1 (data/synthetic/cas1_gatineau_ventes.csv, 15 trimestres 2023T1-2026T3) :
  * l'indice de la médiane et de la moyenne géométrique des prix (biais de composition, décomposition exacte) ;
  * l'indice hédonique à indicatrices temporelles (traite.hedonique.indice_dummies_temps, modèle de référence traite.cas1.matrice
    avec les indicatrices trimestrielles du ch. 18), l'indice chaîné par périodes adjacentes et l'indice d'imputation géométrique
    (égal à l'indice à indicatrices quand les coefficients sont communs) ;
  * un jeu de REVENTES SYNTHÉTIQUES dérivé des 600 ventes (chaque vente observée est la première vente ; la revente est simulée
    avec la loi de génération du jeu : tendance 4,43 %/an, bruit de transaction σ = 0,085 à chaque vente, marche aléatoire
    idiosyncrasique δ² = 0,0004 par mois, 15 % de reventes rapides « rénovées » avec prime) ;
  * les estimateurs de ventes répétées de Bailey-Muth-Nourse (traite.hedonique.ventes_repetees_bmn) et de Case-Shiller pondéré en
    trois étapes (traite.hedonique.ventes_repetees_case_shiller), avec erreurs-types ;
  * l'indice SPAR (prix / valeur au rôle, Bourassa, Hoesli et Sun 2006) construit avec valeur_role_2025 ;
  * la comparaison à la vérité de génération, l'ajustement de conditions du marché du comparable G0404 selon l'indice,
    la sensibilité (millésimes et révision, filtre des reventes rapides, taille de l'échantillon, fenêtre d'estimation) et les
    chiffres des solutions d'exercices.
Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch35_indices.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]))
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parent))
from traite import cas1
from traite import hedonique as h
import ch18_hedonique as c18

ROOT = pathlib.Path(__file__).resolve().parents[3]
GRAINE = 20261031
T = 15                                        # trimestres 2023T1 … 2026T3
LABELS = [f"{2023 + k // 4}T{k % 4 + 1}" for k in range(T)]
BETA_VRAI = c18.VRAIS["mois"]                 # ln(1,0443)/12 : tendance mensuelle de génération
SIGMA = c18.VRAIS["sigma"]                    # bruit de transaction de génération
DELTA2 = 0.0004                               # variance mensuelle de la marche aléatoire idiosyncrasique des reventes (hypothèse)
P_REVENTE = 0.60                              # probabilité qu'une vente observée soit suivie d'une revente dans la fenêtre
PART_FLIP = 0.15                              # part des reventes rapides « rénovées »
PRIME_FLIP = (0.12, 0.04)                     # prime moyenne et écart-type (log) des reventes rapides
FIN_FENETRE = pd.Timestamp("2026-08-31")

# ----------------------------------------------------------------------------- données et vérité
def charger() -> pd.DataFrame:
    v = cas1.charger()["ventes"].copy()
    v["trim"] = c18.trimestre(v["date_vente"])
    v["mois"] = ((v["date_vente"] - cas1.REF_DATE).dt.days / 30.4375).astype(float)
    return v

def mois_mi_trimestre() -> np.ndarray:
    mid = pd.to_datetime([f"{2023 + k // 4}-{3 * (k % 4) + 2:02d}-15" for k in range(T)])
    return ((mid - cas1.REF_DATE).days / 30.4375).to_numpy()

def verite(beta: float = BETA_VRAI) -> np.ndarray:
    """Indice vrai (loi de génération) au milieu de chaque trimestre, base 2023T1 = 100."""
    m = mois_mi_trimestre()
    return 100 * np.exp(beta * (m - m[0]))

# ----------------------------------------------------------------------------- médiane, moyenne, composition
def indices_naifs(v: pd.DataFrame, td: dict) -> pd.DataFrame:
    """Indice de la médiane, de la moyenne arithmétique, de la moyenne géométrique et du prix médian au m² par trimestre ;
    décomposition exacte de la moyenne géométrique : mean(ln P)_t = α + x̄_t'β + τ_t (les résidus MCO somment à zéro par trimestre)."""
    g = v.groupby("trim")
    out = pd.DataFrame({"trimestre": LABELS, "n": g.size().reindex(range(T)).to_numpy(),
                        "mediane": g["prix_vente"].median().to_numpy(), "moyenne": g["prix_vente"].mean().to_numpy(),
                        "moy_geo": np.exp(g["prix_vente"].apply(lambda s: np.log(s).mean())).to_numpy(),
                        "mediane_m2": g.apply(lambda d: (d["prix_vente"] / d["superficie_habitable_m2"]).median()).to_numpy(),
                        "superficie_mediane": g["superficie_habitable_m2"].median().to_numpy(),
                        "part_aylmer": g["secteur"].apply(lambda s: (s == "Aylmer").mean()).to_numpy(),
                        "part_hull": g["secteur"].apply(lambda s: (s == "Hull").mean()).to_numpy(),
                        "age_median": g["age"].median().to_numpy()})
    for c in ("mediane", "moyenne", "moy_geo", "mediane_m2"):
        out[f"indice_{c}"] = 100 * out[c] / out[c].iloc[0]
    # composition : x̄_t'β̂ (hors indicatrices temporelles) relativement à 2023T1
    X = td["X"]; beta = pd.Series(td["mco"]["beta"], index=X.columns)
    cols = [c for c in X.columns if not c.startswith("trim_")]
    xb = (X[cols].values @ beta[cols].values)
    comp = pd.Series(xb, index=v.index).groupby(v["trim"]).mean().to_numpy()
    out["composition_log"] = comp - comp[0]
    out["indice_composition"] = 100 * np.exp(out["composition_log"])
    out["indice_moy_geo_corrigee"] = 100 * np.exp(np.log(out["moy_geo"] / out["moy_geo"].iloc[0]) - out["composition_log"])
    return out

# ----------------------------------------------------------------------------- hédonique : indicatrices, chaînage, imputation
def indice_td(v: pd.DataFrame) -> dict:
    """Indice hédonique à indicatrices trimestrielles (2023T1 = référence) sur le modèle de référence du cas 1."""
    X = c18.matrice_trimestres(v); y = np.log(v["prix_vente"].astype(float)).values
    r = h.mco(X.values, y)
    cols = [c for c in X.columns if c.startswith("trim_")]
    ix = [list(X.columns).index(c) for c in cols]
    tau = r["beta"][ix]; se = r["se_hc1"][ix]
    idx = h.indice_dummies_temps(tau)
    return {"X": X, "y": y, "mco": r, "tau": np.concatenate([[0.0], tau]), "se": np.concatenate([[0.0], se]), "indice": idx,
            "bas": np.concatenate([[100.0], 100 * np.exp(tau - 1.96 * se)]), "haut": np.concatenate([[100.0], 100 * np.exp(tau + 1.96 * se)]),
            "sigma": float(np.sqrt(r["sigma2"])), "r2_ajuste": r["r2_ajuste"], "k": r["k"]}

def indice_chaine(v: pd.DataFrame) -> dict:
    """Indice hédonique chaîné : pour chaque paire de trimestres adjacents (t−1, t), MCO sur les ventes des deux trimestres avec une
    indicatrice de t (coefficients des caractéristiques propres à la paire) ; I_t = I_{t−1} exp(δ_t)."""
    delta = np.zeros(T); se = np.zeros(T)
    for t in range(1, T):
        sub = v[v["trim"].isin([t - 1, t])]
        X = cas1.matrice(sub).drop(columns=["mois"])
        X = X[[c for c in X.columns if X[c].std() > 0 or c == "const"]]
        X["d_t"] = (sub["trim"] == t).astype(float)
        r = h.mco(X.values, np.log(sub["prix_vente"].astype(float)).values)
        j = list(X.columns).index("d_t"); delta[t] = r["beta"][j]; se[t] = r["se_hc1"][j]
    ln_idx = np.cumsum(delta)
    return {"delta": delta, "se_delta": se, "indice": 100 * np.exp(ln_idx), "se_cumul": np.sqrt(np.cumsum(se ** 2))}

def indice_imputation(v: pd.DataFrame, td: dict) -> dict:
    """Imputation géométrique (Laspeyres : panier de la période de base ; Paasche : panier courant ; Fisher = moyenne géométrique)
    avec le modèle à indicatrices (coefficients communs) : chaque vente de la période 0 (resp. t) est « revendue » à la période t (resp. 0)
    au prix prédit. Avec des coefficients communs, les trois indices coïncident avec exp(τ_t) ; la fonction le vérifie numériquement."""
    X = td["X"]; beta = td["mco"]["beta"]; cols = list(X.columns)
    tcols = [c for c in cols if c.startswith("trim_")]
    def predire_a(Xsub: pd.DataFrame, t: int) -> np.ndarray:
        Z = Xsub.copy()
        for c in tcols: Z[c] = 0.0
        if t > 0: Z[f"trim_{t}"] = 1.0
        return Z.values @ beta
    lasp = np.zeros(T); paas = np.zeros(T)
    X0 = X[v["trim"] == 0]
    for t in range(T):
        Xt = X[v["trim"] == t]
        lasp[t] = np.mean(predire_a(X0, t) - predire_a(X0, 0))
        paas[t] = np.mean(predire_a(Xt, t) - predire_a(Xt, 0))
    fisher = 0.5 * (lasp + paas)
    return {"laspeyres": 100 * np.exp(lasp), "paasche": 100 * np.exp(paas), "fisher": 100 * np.exp(fisher),
            "ecart_max_vs_td": float(np.max(np.abs(100 * np.exp(fisher) - td["indice"])))}

# ----------------------------------------------------------------------------- reventes synthétiques
def generer_reventes(v: pd.DataFrame, graine: int = GRAINE, p_revente: float = P_REVENTE, part_flip: float = PART_FLIP,
                     delta2: float = DELTA2, beta: float = BETA_VRAI, sigma: float = SIGMA) -> pd.DataFrame:
    """Paires de reventes synthétiques. Pour chaque vente observée (première vente, prix P1 à la date t1), avec probabilité p_revente :
    intervalle Δ (mois) tiré d'une loi déclarée (reventes « rapides » : 3 à 9 mois, prime log N(0,08 ; 0,03) ; autres : 6 + Exp(16), ≥ 6) ;
    ln P2 = ln P1 + β Δ + prime·flip + w + ε2 − ε1', avec ε1', ε2 ~ N(0, σ²) (bruit de transaction à chaque vente) et w ~ N(0, δ² Δ).
    Seules les reventes datées au plus tard le 31 août 2026 sont conservées ; prix arrondis à 500 $."""
    rng = np.random.default_rng(graine)
    rows = []
    for _, s in v.iterrows():
        if rng.random() > p_revente:
            continue
        flip = rng.random() < part_flip
        if flip:
            gap = int(rng.integers(3, 10)); prime = rng.normal(*PRIME_FLIP)
        else:
            gap = int(max(6, round(6 + rng.exponential(16.0)))); prime = 0.0
        d2 = s["date_vente"] + pd.DateOffset(months=gap)
        if d2 > FIN_FENETRE:
            continue
        w = rng.normal(0.0, math.sqrt(delta2 * gap)); e1 = rng.normal(0.0, sigma); e2 = rng.normal(0.0, sigma)
        lnp2 = math.log(s["prix_vente"]) + beta * gap + prime + w + e2 - e1
        p2 = 500 * round(math.exp(lnp2) / 500)
        rows.append({"id_vente": s["id_vente"], "secteur": s["secteur"], "date_1": s["date_vente"], "date_2": d2, "prix_1": float(s["prix_vente"]),
                     "prix_2": float(p2), "t_1": int(s["trim"]), "t_2": int(c18.trimestre(pd.Series([d2]))[0]), "intervalle_mois": gap,
                     "flip": bool(flip), "prime": prime})
    df = pd.DataFrame(rows)
    df["ln_ratio"] = np.log(df["prix_2"] / df["prix_1"])
    df["id_paire"] = [f"R{i + 1:04d}" for i in range(len(df))]
    return df

def resume_reventes(p: pd.DataFrame) -> dict:
    return {"n_paires": len(p), "n_flips": int(p["flip"].sum()), "intervalle_median": float(p["intervalle_mois"].median()),
            "intervalle_moyen": float(p["intervalle_mois"].mean()), "intervalle_min": int(p["intervalle_mois"].min()), "intervalle_max": int(p["intervalle_mois"].max()),
            "part_moins_12": float((p["intervalle_mois"] < 12).mean()), "variation_mediane": float(np.exp(p["ln_ratio"].median()) - 1),
            "variation_moyenne_log": float(p["ln_ratio"].mean()), "sd_ln_ratio": float(p["ln_ratio"].std(ddof=1)),
            "paires_par_annee_1": p["date_1"].dt.year.value_counts().sort_index().to_dict(),
            "paires_par_annee_2": p["date_2"].dt.year.value_counts().sort_index().to_dict(),
            "meme_trimestre": int((p["t_1"] == p["t_2"]).sum())}

# ----------------------------------------------------------------------------- ventes répétées : BMN et Case–Shiller avec erreurs-types
def _se_mco_sans_constante(D: np.ndarray, y: np.ndarray, beta: np.ndarray, w: np.ndarray | None = None) -> np.ndarray:
    if w is None: w = np.ones(len(y))
    Dw = D * w[:, None]; yw = y * w
    res = yw - Dw @ beta; n, k = D.shape
    s2 = res @ res / (n - k)
    return np.sqrt(np.diag(s2 * np.linalg.pinv(Dw.T @ Dw)))

def bmn(p: pd.DataFrame) -> dict:
    r = h.ventes_repetees_bmn(p["prix_1"], p["prix_2"], p["t_1"], p["t_2"], T)
    se = _se_mco_sans_constante(r["D"], r["y"], r["beta"])
    r["se"] = np.concatenate([[0.0], se]); r["bas"] = 100 * np.exp(np.concatenate([[0.0], r["beta"]]) - 1.96 * r["se"])
    r["haut"] = 100 * np.exp(np.concatenate([[0.0], r["beta"]]) + 1.96 * r["se"]); r["sigma"] = float(np.sqrt(r["residus"] @ r["residus"] / (len(r["y"]) - T + 1)))
    # vérification de la forme fermée : β = (D'D)^{-1} D'y
    r["beta_ferme"] = np.linalg.solve(r["D"].T @ r["D"], r["D"].T @ r["y"])
    return r

def case_shiller(p: pd.DataFrame) -> dict:
    r = h.ventes_repetees_case_shiller(p["prix_1"], p["prix_2"], p["t_1"], p["t_2"], T)
    e1 = h.ventes_repetees_bmn(p["prix_1"], p["prix_2"], p["t_1"], p["t_2"], T)
    gap_t = (p["t_2"] - p["t_1"]).to_numpy(float)
    se = _se_mco_sans_constante(e1["D"], e1["y"], r["beta"], r["poids"])
    r["se"] = np.concatenate([[0.0], se]); r["bas"] = 100 * np.exp(np.concatenate([[0.0], r["beta"]]) - 1.96 * r["se"])
    r["haut"] = 100 * np.exp(np.concatenate([[0.0], r["beta"]]) + 1.96 * r["se"])
    r["gap_trimestres"] = gap_t; r["residus_bmn2"] = e1["residus"] ** 2
    a = r["coef_variance"]                      # var(ε) = a0 + a1·(intervalle en trimestres)
    r["variance_a_gap"] = {g: float(a[0] + a[1] * g) for g in (1, 2, 4, 8, 12)}
    r["poids_relatif_a_gap"] = {g: float(math.sqrt((a[0] + a[1] * 1) / (a[0] + a[1] * g))) for g in (1, 2, 4, 8, 12)}
    # pente théorique : var = 2σ² + δ²·(3 mois par trimestre)·gap, en unités de trimestre
    r["variance_theorique"] = (2 * SIGMA ** 2, DELTA2 * 3)
    return r

# ----------------------------------------------------------------------------- SPAR
def spar(v: pd.DataFrame) -> dict:
    """Indice SPAR (sale price appraisal ratio) : arithmétique [(ΣP_t/ΣA_t)/(ΣP_0/ΣA_0)] et géométrique [exp(mean ln(P/A)_t − mean ln(P/A)_0)],
    A = valeur au rôle 2025 (date de référence 2024-07-01, fictive)."""
    g = v.groupby("trim")
    ratio_arith = (g["prix_vente"].sum() / g["valeur_role_2025"].sum()).to_numpy()
    lr = np.log(v["prix_vente"] / v["valeur_role_2025"])
    ratio_geo = lr.groupby(v["trim"]).mean().to_numpy(); sd_geo = lr.groupby(v["trim"]).std(ddof=1).to_numpy(); n = g.size().to_numpy()
    se_geo = np.sqrt(sd_geo ** 2 / n + sd_geo[0] ** 2 / n[0]); se_geo[0] = 0.0
    return {"ratio_arith": ratio_arith, "ratio_geo": np.exp(ratio_geo), "indice_arith": 100 * ratio_arith / ratio_arith[0],
            "indice_geo": 100 * np.exp(ratio_geo - ratio_geo[0]), "se_log": se_geo,
            "bas": 100 * np.exp(ratio_geo - ratio_geo[0] - 1.96 * se_geo), "haut": 100 * np.exp(ratio_geo - ratio_geo[0] + 1.96 * se_geo),
            "sd_ln_ratio": float(lr.std(ddof=1)), "ratio_median_global": float(np.median(v["prix_vente"] / v["valeur_role_2025"]))}

# ----------------------------------------------------------------------------- comparaison
def precision(idx: np.ndarray, vrai: np.ndarray) -> dict:
    e = 100 * (idx / vrai - 1)
    return {"EAM_pct": float(np.mean(np.abs(e[1:]))), "RMSE_pct": float(np.sqrt(np.mean(e[1:] ** 2))), "max_abs_pct": float(np.max(np.abs(e))),
            "croissance_totale_pct": float(idx[-1] - 100), "ecart_final_pct": float(e[-1]),
            "variations_trim_sd_pct": float(np.std(np.diff(np.log(idx))) * 100)}

def tableau_indices(vrai, naifs, td, ch, b, cs, sp) -> pd.DataFrame:
    return pd.DataFrame({"trimestre": LABELS, "n_ventes": naifs["n"], "vrai": vrai, "mediane": naifs["indice_mediane"].to_numpy(),
                         "moy_geo": naifs["indice_moy_geo"].to_numpy(), "hedonique": td["indice"], "hed_bas": td["bas"], "hed_haut": td["haut"],
                         "chaine": ch["indice"], "bmn": b["indice"], "bmn_bas": b["bas"], "bmn_haut": b["haut"], "cs": cs["indice"],
                         "cs_bas": cs["bas"], "cs_haut": cs["haut"], "spar": sp["indice_geo"], "spar_arith": sp["indice_arith"]})

def ajustement_comparable(v: pd.DataFrame, tab: pd.DataFrame, ref: dict, id_vente: str = "G0404", date_eval: str = "2026-09-01") -> dict:
    """Ajustement de conditions du marché d'un comparable du ch. 16 selon l'indice retenu : A = P × (I_eval / I_vente − 1)."""
    c = v[v["id_vente"] == id_vente].iloc[0]; P = float(c["prix_vente"]); t1 = int(c["trim"]); t2 = T - 1
    dm = (pd.Timestamp(date_eval) - c["date_vente"]).days / 30.4375
    bm = ref["beta"][list(ref["X"].columns).index("mois")]
    out = {"prix": P, "trimestre_vente": LABELS[t1], "mois_ecoules": float(dm), "tendance": P * (math.exp(bm * dm) - 1)}
    for k in ("mediane", "hedonique", "chaine", "bmn", "cs", "spar"):
        out[k] = P * (tab[k].iloc[t2] / tab[k].iloc[t1] - 1)
    out["vrai"] = P * (math.exp(BETA_VRAI * dm) - 1)
    return out

# ----------------------------------------------------------------------------- sensibilité
def sensibilite(v: pd.DataFrame, p: pd.DataFrame, vrai: np.ndarray, td: dict, graine: int = GRAINE) -> dict:
    out = {}
    # (a) millésimes : révision des indices de ventes répétées quand de nouvelles paires s'ajoutent
    mill = {}
    for fin, lab in ((9, "2025T2"), (11, "2025T4"), (T - 1, "2026T3")):
        sub = p[p["t_2"] <= fin]
        rb = h.ventes_repetees_bmn(sub["prix_1"], sub["prix_2"], sub["t_1"], sub["t_2"], fin + 1)
        rc = h.ventes_repetees_case_shiller(sub["prix_1"], sub["prix_2"], sub["t_1"], sub["t_2"], fin + 1)
        mill[lab] = {"n": len(sub), "bmn": rb["indice"], "cs": rc["indice"]}
    out["millesimes"] = mill
    rev = {}
    for t, lab in ((9, "2025T2"), (7, "2024T4"), (5, "2024T2")):
        rev[lab] = {"bmn_v2025T2": float(mill["2025T2"]["bmn"][t]), "bmn_v2025T4": float(mill["2025T4"]["bmn"][t]), "bmn_v2026T3": float(mill["2026T3"]["bmn"][t]),
                    "cs_v2025T2": float(mill["2025T2"]["cs"][t]), "cs_v2026T3": float(mill["2026T3"]["cs"][t])}
    out["revisions"] = rev
    # (b) filtre des reventes rapides
    filt = {}
    for seuil, lab in ((0, "toutes"), (12, "intervalle ≥ 12 mois"), (18, "intervalle ≥ 18 mois")):
        sub = p[p["intervalle_mois"] >= seuil] if seuil else p
        rb = h.ventes_repetees_bmn(sub["prix_1"], sub["prix_2"], sub["t_1"], sub["t_2"], T)
        rc = h.ventes_repetees_case_shiller(sub["prix_1"], sub["prix_2"], sub["t_1"], sub["t_2"], T)
        filt[lab] = {"n": len(sub), "n_flips": int(sub["flip"].sum()), "bmn": precision(rb["indice"], vrai), "cs": precision(rc["indice"], vrai),
                     "bmn_indice": rb["indice"], "cs_indice": rc["indice"]}
    sans_flip = p[~p["flip"]]
    rb = h.ventes_repetees_bmn(sans_flip["prix_1"], sans_flip["prix_2"], sans_flip["t_1"], sans_flip["t_2"], T)
    filt["sans les reventes rénovées (oracle)"] = {"n": len(sans_flip), "n_flips": 0, "bmn": precision(rb["indice"], vrai), "bmn_indice": rb["indice"]}
    out["filtre"] = filt
    # (c) taille de l'échantillon : sous-échantillons aléatoires de paires, erreur-type du dernier trimestre et EAM (moyenne de 200 tirages)
    rng = np.random.default_rng(graine); taille = {}
    for frac in (1.0, 0.5, 0.25):
        eam, fin = [], []
        for _ in range(200 if frac < 1 else 1):
            sub = p.sample(frac=frac, random_state=int(rng.integers(1e9))) if frac < 1 else p
            try:
                rb = h.ventes_repetees_bmn(sub["prix_1"], sub["prix_2"], sub["t_1"], sub["t_2"], T)
            except Exception:
                continue
            eam.append(precision(rb["indice"], vrai)["EAM_pct"]); fin.append(rb["indice"][-1])
        taille[frac] = {"n": int(round(frac * len(p))), "EAM_moyen_pct": float(np.mean(eam)), "sd_indice_final": float(np.std(fin)) if len(fin) > 1 else float("nan")}
    out["taille"] = taille
    # (d) fenêtre d'estimation de l'indice hédonique : huit derniers trimestres seulement, rechaîné sur 2024T4
    sub = v[v["trim"] >= 7]
    X = cas1.matrice(sub).drop(columns=["mois"]); X = X[[c for c in X.columns if X[c].std() > 0 or c == "const"]]
    for k in range(8, T): X[f"trim_{k}"] = (sub["trim"] == k).astype(float)
    r = h.mco(X.values, np.log(sub["prix_vente"].astype(float)).values)
    tau = np.array([r["beta"][list(X.columns).index(f"trim_{k}")] for k in range(8, T)])
    idx_court = td["indice"][7] * np.exp(np.concatenate([[0.0], tau]))
    out["fenetre"] = {"n": len(sub), "indice_court": idx_court, "indice_long": td["indice"][7:], "ecart_max_pct": float(np.max(np.abs(idx_court / td["indice"][7:] - 1)) * 100),
                      "sigma_court": float(np.sqrt(r["sigma2"]))}
    return out

# ----------------------------------------------------------------------------- exercices
def exercices(tab: pd.DataFrame, b: dict, cs: dict, sp: dict, ref: dict) -> dict:
    e = {}
    # 35.1 : médiane et composition (deux trimestres fictifs)
    e["1"] = {"mediane_t0": 450_000, "mediane_t1": 480_000, "var_mediane_pct": (480_000 / 450_000 - 1) * 100,
              "composition_pct": (math.exp(0.539 * math.log(150 / 140)) - 1) * 100,
              "variation_pure_pct": (480_000 / 450_000 / math.exp(0.539 * math.log(150 / 140)) - 1) * 100}
    # 35.2 : indice à indicatrices
    tau = np.array([0.0, 0.031, 0.058, 0.102]); e["2"] = {"indice": 100 * np.exp(tau), "variation_t3_t2_pct": (math.exp(0.102 - 0.058) - 1) * 100,
                                                          "ic_t3": (100 * math.exp(0.102 - 1.96 * 0.021), 100 * math.exp(0.102 + 1.96 * 0.021))}
    # 35.3 : BMN à la main sur quatre paires et trois périodes
    D = np.array([[1.0, 0.0], [-1.0, 1.0], [0.0, 1.0], [1.0, 0.0]]); y = np.log(np.array([462_000 / 430_000, 512_000 / 470_000, 540_000 / 495_000, 455_000 / 436_000]))
    beta = np.linalg.solve(D.T @ D, D.T @ y)
    e["3"] = {"y": y, "DtD": D.T @ D, "Dty": D.T @ y, "beta": beta, "indice": 100 * np.exp(np.concatenate([[0.0], beta])), "residus": y - D @ beta}
    # 35.4 : Case–Shiller : variance et poids
    a0, a1 = 0.0145, 0.0012; gaps = np.array([1, 4, 8, 12])
    var = a0 + a1 * gaps; e["4"] = {"var": var, "poids": 1 / np.sqrt(var), "poids_rel": np.sqrt(var[0] / var), "sd_pct": 100 * np.sqrt(var)}
    # 35.5 : SPAR sur dix ventes fictives
    P = np.array([512_000, 468_500, 601_000, 445_000, 530_000, 487_500, 575_000, 455_000, 498_000, 522_000.0]); A = np.array([470_000, 452_000, 545_000, 425_000, 480_000, 468_000, 510_000, 440_000, 455_000, 471_000.0])
    base_ratio = 1.045; r_ar = P.sum() / A.sum(); r_geo = math.exp(np.mean(np.log(P / A)))
    e["5"] = {"ratio_arith": r_ar, "ratio_geo": r_geo, "spar_arith": 100 * r_ar / base_ratio, "spar_geo": 100 * r_geo / base_ratio, "ratios": P / A}
    # 35.6 : ajustement de G0404 selon l'indice
    e["6"] = {"P": 506_500.0, "I_2025T2": {k: float(tab[k].iloc[9]) for k in ("hedonique", "bmn", "cs", "mediane")},
              "I_2026T3": {k: float(tab[k].iloc[14]) for k in ("hedonique", "bmn", "cs", "mediane")}}
    e["6"]["ajustements"] = {k: 506_500.0 * (e["6"]["I_2026T3"][k] / e["6"]["I_2025T2"][k] - 1) for k in e["6"]["I_2025T2"]}
    bm = ref["beta"][list(ref["X"].columns).index("mois")]; e["6"]["tendance"] = 506_500.0 * (math.exp(bm * 15.2) - 1)
    # 35.7 : révision et biais de sélection (chiffres du script : voir sensibilité) — valeurs conceptuelles
    # 35.8 : script (pas de chiffre fixe)
    return e

# ----------------------------------------------------------------------------- principal
def main(verbose: bool = True) -> dict:
    v = charger(); ref = cas1.hedonique(v.drop(columns=["trim", "mois"]))
    vrai = verite()
    td = indice_td(v); ch = indice_chaine(v); imp = indice_imputation(v, td)
    naifs = indices_naifs(v, td)
    p = generer_reventes(v); rp = resume_reventes(p)
    b = bmn(p); cs = case_shiller(p); sp = spar(v)
    tab = tableau_indices(vrai, naifs, td, ch, b, cs, sp)
    prec = {k: precision(tab[k].to_numpy(), vrai) for k in ("mediane", "moy_geo", "hedonique", "chaine", "bmn", "cs", "spar", "spar_arith")}
    prec["hedonique"]["largeur_ic_finale_pct"] = float((td["haut"][-1] / td["bas"][-1] - 1) * 100)
    prec["bmn"]["largeur_ic_finale_pct"] = float((b["haut"][-1] / b["bas"][-1] - 1) * 100)
    prec["cs"]["largeur_ic_finale_pct"] = float((cs["haut"][-1] / cs["bas"][-1] - 1) * 100)
    prec["spar"]["largeur_ic_finale_pct"] = float((sp["haut"][-1] / sp["bas"][-1] - 1) * 100)
    # corrélation entre l'écart médiane − hédonique et la composition
    ecart_med = np.log(naifs["indice_mediane"].to_numpy() / td["indice"]); corr_comp = float(np.corrcoef(ecart_med, naifs["composition_log"])[0, 1])
    corr_surf = float(np.corrcoef(ecart_med, np.log(naifs["superficie_mediane"]))[0, 1])
    adj = ajustement_comparable(v, tab, ref)
    sens = sensibilite(v, p, vrai, td)
    ex = exercices(tab, b, cs, sp, ref)
    out = {"ventes": v, "ref": ref, "vrai": vrai, "td": td, "chaine": ch, "imputation": imp, "naifs": naifs, "paires": p, "resume_paires": rp,
           "bmn": b, "cs": cs, "spar": sp, "tableau": tab, "precision": prec, "corr_ecart_composition": corr_comp, "corr_ecart_superficie": corr_surf,
           "ajustement": adj, "sensibilite": sens, "exercices": ex, "beta_vrai": BETA_VRAI}
    if verbose:
        pd.set_option("display.width", 240); pd.set_option("display.max_columns", 40)
        print(f"== Vérité de génération : β = {BETA_VRAI:.6f}/mois ({(math.exp(12 * BETA_VRAI) - 1) * 100:.2f} %/an) ; indice final {vrai[-1]:.2f}")
        print("== Médiane, moyennes et composition :\n", naifs.round(3).to_string(index=False))
        print(f"   corrélation(ln médiane/hédonique, composition) = {corr_comp:.3f} ; avec ln superficie médiane = {corr_surf:.3f}")
        print(f"== Hédonique à indicatrices : k = {td['k']}, σ = {td['sigma']:.4f}, R² aj. = {td['r2_ajuste']:.4f}")
        print(f"== Imputation géométrique : écart max à l'indice à indicatrices = {imp['ecart_max_vs_td']:.2e}")
        print("== Reventes synthétiques :", {k: (round(x, 4) if isinstance(x, float) else x) for k, x in rp.items()})
        print(f"== BMN : σ = {b['sigma']:.4f} ; écart max forme fermée vs lstsq = {np.max(np.abs(b['beta'] - b['beta_ferme'])):.1e}")
        print(f"== Case–Shiller : var(ε) = {cs['coef_variance'][0]:.5f} + {cs['coef_variance'][1]:.5f} × intervalle (trimestres) ; théorie {cs['variance_theorique']}")
        print("   poids relatifs :", {k: round(x, 3) for k, x in cs["poids_relatif_a_gap"].items()})
        print(f"== SPAR : sd ln(P/A) = {sp['sd_ln_ratio']:.4f} ; ratio médian global {sp['ratio_median_global']:.4f}")
        print("== Indices par trimestre :\n", tab.round(2).to_string(index=False))
        print("== Précision vs vérité :")
        for k, d in prec.items(): print(f"   {k:11s}", {kk: round(x, 2) for kk, x in d.items()})
        print("== Ajustement de G0404 :", {k: (round(x, 1) if isinstance(x, float) else x) for k, x in adj.items()})
        print("== Sensibilité — révisions :", sens["revisions"])
        for lab, d in sens["filtre"].items(): print(f"   filtre {lab:36s} n={d['n']:3d} flips={d['n_flips']:2d} BMN EAM={d['bmn']['EAM_pct']:.2f} final={d['bmn']['croissance_totale_pct']:.2f}", f"CS EAM={d['cs']['EAM_pct']:.2f}" if "cs" in d else "")
        print("   taille :", sens["taille"])
        print("   fenêtre hédonique :", {k: (np.round(x, 2) if isinstance(x, np.ndarray) else round(x, 3)) for k, x in sens["fenetre"].items()})
        print("== Exercices :")
        for k, e in ex.items(): print("  ", k, {kk: (np.round(vv, 4) if isinstance(vv, (np.ndarray, float)) else vv) for kk, vv in e.items()})
    return out

if __name__ == "__main__":
    main()