#!/usr/bin/env python3
"""Chapitre 5 — Économie urbaine et rente foncière : calculs de l'exemple intégral et des encadrés.

Imprime tous les chiffres cités dans le chapitre 5 :
  * modèle monocentrique (Alonso–Muth–Mills) simulé avec des PARAMÈTRES PÉDAGOGIQUES DÉCLARÉS (section PARAMS) :
    gradient du prix du logement p(d) = p0 (1 - t d / y)^(1/beta), gradient exponentiel approché, rente foncière
    r(d) = r0 (p(d)/p0)^(1/alpha), densité, frontière urbaine, courbes d'enchère de von Thünen et d'Alonso ;
  * gradient prix/distance observé dans le dossier du cas 1 (Gatineau, DONNÉES SYNTHÉTIQUES, graine 20261031) :
    tranches de distance, régressions avec et sans indicatrices de secteur (traite.hedonique.mco), coefficient
    `dist_centre` du modèle hédonique de référence (traite.cas1.hedonique) ;
  * décomposition interrégionale de la valeur par logement du cas 4 (traite.cas4.sommaire) ;
  * valeur résiduelle du terrain et effet de levier ;
  * capitalisation de la taxe foncière : formule exacte et approchée, taux effectif, taux 2026 de Québec,
    Gatineau et Montréal (sources : notes-B/E du dépôt), grille de sensibilité.

Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch05_rente.py
Les régressions utilisent traite.hedonique.mco ; les données des cas passent par traite.cas1 et traite.cas4."""
from __future__ import annotations
import sys, pathlib, math
import numpy as np, pandas as pd
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
from traite import cas1, cas4
from traite import hedonique as h

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

# ----------------------------------------------------------------------------- paramètres pédagogiques (déclarés)
PARAMS = {
    "y": 80000.0,     # revenu annuel du ménage ($/an)
    "t": 700.0,       # coût annuel de navettage par kilomètre de distance au centre ($/km/an), temps compris
    "beta": 0.30,     # part du revenu net de transport consacrée au logement
    "alpha": 0.30,    # part du terrain dans le coût de production du logement (Cobb-Douglas)
    "p0": 260.0,      # prix annuel du service de logement au centre ($/m² de plancher/an)
    "r0": 80.0,       # rente foncière annuelle au centre ($/m² de terrain/an)
    "ra": 2.0,        # rente alternative (agricole, conversion comprise) à la frontière ($/m²/an)
    "i": 0.05,        # taux de capitalisation des rentes foncières (annuel)
}
# courbes d'enchère linéaires de trois usages (illustration d'Alonso) : (rente au centre $/m²/an, pente $/m²/an/km)
ALONSO = {"Commerces et bureaux": (300.0, 25.0), "Résidentiel dense": (120.0, 6.0), "Résidentiel pavillonnaire": (60.0, 1.8)}
# von Thünen (exercice 5.1) : rendement Y (t/ha), prix p ($/t), coût c ($/t), fret f ($/t/km)
THUNEN = {"Maraîchage": dict(Y=20.0, p=300.0, c=250.0, f=3.0), "Grandes cultures": dict(Y=4.0, p=500.0, c=350.0, f=5.0)}

# taux de taxe foncière 2026 (par 100 $ d'évaluation) — sources : vqtaux2026 ; gatineaureglement994 ; ledevoirmtl2026 (secondaire) ; taxescolaire2026
TAUX_2026 = {
    "Québec (1 à 5 logements)": 0.7464, "Québec (6 logements et plus)": 0.7697,
    "Gatineau (taux de base)": 0.5224, "Montréal (taxe foncière générale, hors arrondissement)": 0.4631,
    "Taxe scolaire 2026-2027 (toute la province)": 0.07899,
}
FACTEUR_2026 = {"Gatineau": 1.02, "Montréal": 1.00}   # facteurs comparatifs 2026 publiés (gatmutations ; mtlmutations2026)


# ----------------------------------------------------------------------------- von Thünen et Alonso
def rente_thunen(d, Y: float, p: float, c: float, f: float) -> np.ndarray:
    """R(d) = Y (p - c) - Y f d  (rente par hectare et par an)."""
    return Y * (p - c) - Y * f * np.asarray(d, float)

def distance_limite_thunen(Y: float, p: float, c: float, f: float) -> float:
    """Distance où la rente s'annule : d* = (p - c) / f."""
    return (p - c) / f

def frontiere_usages(R1: float, t1: float, R2: float, t2: float) -> float:
    """Distance où deux courbes d'enchère linéaires R1 - t1 d et R2 - t2 d se croisent."""
    return (R1 - R2) / (t1 - t2)

def enveloppe_alonso(d, usages: dict = ALONSO, ra: float = PARAMS["ra"]) -> tuple[np.ndarray, list]:
    d = np.asarray(d, float)
    courbes = {k: R0 - t * d for k, (R0, t) in usages.items()}
    env = np.maximum.reduce([*courbes.values(), np.full_like(d, ra)])
    return env, courbes


# ----------------------------------------------------------------------------- modèle monocentrique
def prix_logement(d, P: dict = PARAMS) -> np.ndarray:
    """p(d) = p0 (1 - t d / y)^(1/beta) : solution exacte de la condition de Muth avec part budgétaire constante."""
    d = np.asarray(d, float)
    base = np.clip(1.0 - P["t"] * d / P["y"], 1e-12, None)
    return P["p0"] * base ** (1.0 / P["beta"])

def gradient_exponentiel(P: dict = PARAMS) -> float:
    """Pente logarithmique du prix au centre : gamma = t / (beta y) (approximation exponentielle)."""
    return P["t"] / (P["beta"] * P["y"])

def prix_exponentiel(d, P: dict = PARAMS) -> np.ndarray:
    return P["p0"] * np.exp(-gradient_exponentiel(P) * np.asarray(d, float))

def consommation_logement(d, P: dict = PARAMS) -> np.ndarray:
    """q(d) = beta (y - t d) / p(d) : surface de plancher consommée par ménage (m²)."""
    d = np.asarray(d, float)
    return P["beta"] * (P["y"] - P["t"] * d) / prix_logement(d, P)

def rente_fonciere(d, P: dict = PARAMS) -> np.ndarray:
    """r(d) = r0 (p(d)/p0)^(1/alpha) : profit nul des producteurs de logement (Cobb-Douglas)."""
    return P["r0"] * (prix_logement(d, P) / P["p0"]) ** (1.0 / P["alpha"])

def intensite_capital(d, P: dict = PARAMS) -> np.ndarray:
    """Plancher par m² de terrain, normalisé à 1 au centre : h(d)/h(0) = (p(d)/p0)^((1-alpha)/alpha)."""
    return (prix_logement(d, P) / P["p0"]) ** ((1.0 - P["alpha"]) / P["alpha"])

def densite(d, P: dict = PARAMS) -> np.ndarray:
    """Ménages par m² de terrain, normalisée à 1 au centre : D(d) = h(d)/q(d)."""
    d = np.asarray(d, float)
    D = intensite_capital(d, P) / consommation_logement(d, P)
    return D / (intensite_capital(0.0, P) / consommation_logement(0.0, P))

def frontiere_urbaine(P: dict = PARAMS) -> float:
    """Distance d_bar où r(d_bar) = ra (solution analytique)."""
    ratio_p = (P["ra"] / P["r0"]) ** P["alpha"]            # p(d_bar)/p0
    return (1.0 - ratio_p ** P["beta"]) * P["y"] / P["t"]

def demi_distance(gamma: float) -> float:
    """Distance qui divise le prix par deux sous un gradient exponentiel : ln 2 / gamma."""
    return math.log(2.0) / gamma

def valeur_terrain(rente: float, i: float = PARAMS["i"], g: float = 0.0) -> float:
    """Valeur du terrain = rente / (i - g) (rente perpétuelle, croissance g)."""
    return rente / (i - g)

def decomposition_capozza(ra: float, i: float, cout_conversion: float, rente_d: float, prime_croissance: float) -> dict:
    """Capozza–Helsley (1989) : P_L = ra/i + C + (r(d) - ra)/i + prime de croissance."""
    comp = {"valeur_agricole": ra / i, "cout_conversion": cout_conversion, "accessibilite": (rente_d - ra) / i, "prime_croissance": prime_croissance}
    comp["total"] = sum(comp.values()); return comp


# ----------------------------------------------------------------------------- rente résiduelle et levier
def valeur_residuelle(V: float, CN: float, depreciation: float) -> float:
    """V_L = V - (C_N - D) : valeur du terrain par résidu."""
    return V - (CN - depreciation)

def levier_residuel(V: float, VL: float) -> float:
    """Élasticité de V_L par rapport à V : V / V_L."""
    return V / VL

def elasticite_rente_prix(alpha: float) -> float:
    return 1.0 / alpha


# ----------------------------------------------------------------------------- capitalisation fiscale
def valeur_avec_taxe(R: float, i: float, tau: float) -> float:
    """V = R / (i + tau) : rente nette de taxe capitalisée (taxe proportionnelle à la valeur)."""
    return R / (i + tau)

def capitalisation_taxe(V0: float, tau0: float, tau1: float, i: float, kappa: float = 1.0) -> dict:
    """Effet d'un changement de taux de taxe tau0 -> tau1 sur une valeur V0 établie sous tau0.
    exact : V1 = R/(i+tau1) avec R = V0 (i+tau0) ; approché : dV = -kappa * dT / i avec dT = (tau1-tau0) V0."""
    R = V0 * (i + tau0)
    V1 = R / (i + tau1)
    dT = (tau1 - tau0) * V0
    return {"R": R, "V1_exact": V1, "dV_exact": V1 - V0, "dT": dT, "dV_approx": -kappa * dT / i,
            "dV_exact_kappa": kappa * (V1 - V0), "dV_pct_exact": V1 / V0 - 1}

def taux_effectif(tau_role: float, facteur: float) -> float:
    """Taux sur la valeur marchande : tau_eff = tau_role / F, où F est le facteur comparatif (V_marché = F × V_rôle)."""
    return tau_role / facteur

def grille_capitalisation(V0: float, d_taux_pb: list, taux_cap: list, kappa: float = 1.0) -> pd.DataFrame:
    """Perte de valeur (en $) pour un écart de taux de taxe (points de base de la valeur) selon le taux de capitalisation."""
    rows = []
    for pb in d_taux_pb:
        row = {"ecart_pb": pb}
        for r in taux_cap:
            row[f"r={r:.0%}"] = -kappa * (pb / 10000.0) * V0 / r
        rows.append(row)
    return pd.DataFrame(rows).set_index("ecart_pb")


# ----------------------------------------------------------------------------- cas 1 : gradient observé
TRANCHES = [0, 2, 4, 6, 8, 10, 12, 20, 35]
ETIQUETTES = ["0-2", "2-4", "4-6", "6-8", "8-10", "10-12", "12-20", "20 et plus"]

def tranches_distance(v: pd.DataFrame) -> pd.DataFrame:
    v = v.copy(); v["prix_m2"] = v["prix_vente"] / v["superficie_habitable_m2"]
    v["tranche"] = pd.cut(v["dist_centre_km"], TRANCHES, labels=ETIQUETTES, right=True)
    g = v.groupby("tranche", observed=True).agg(n=("prix_vente", "size"), dist_mediane=("dist_centre_km", "median"),
                                                 prix_median=("prix_vente", "median"), prix_m2_median=("prix_m2", "median"),
                                                 superficie_mediane=("superficie_habitable_m2", "median"),
                                                 terrain_median=("superficie_terrain_m2", "median"))
    g["secteur_dominant"] = v.groupby("tranche", observed=True)["secteur"].agg(lambda s: s.value_counts().index[0])
    return g

def gradients_cas1(v: pd.DataFrame, modele: dict | None = None) -> dict:
    """Quatre spécifications du gradient prix/distance (coefficients en log par km)."""
    y = np.log(v["prix_vente"].astype(float).values)
    d = v["dist_centre_km"].astype(float).values
    out = {}
    # (a) brut : ln(P/S) = a - gamma d
    r = h.mco(np.column_stack([np.ones_like(d), d]), np.log(v["prix_vente"] / v["superficie_habitable_m2"]).values)
    out["brut_prix_m2"] = {"gamma": -r["beta"][1], "se": r["se_hc1"][1], "r2": r["r2"], "n": r["n"]}
    # (b) contrôlé : caractéristiques du modèle de référence sans secteurs, distance linéaire
    X = cas1.matrice(v).drop(columns=["dist_centre"] + [f"sect_{s}" for s in cas1.SECTEURS[1:]])
    X["dist"] = d
    r = h.mco(X.values, y)
    out["controle_lineaire"] = {"gamma": -r["beta"][-1], "se": r["se_hc1"][-1], "r2": r["r2"], "n": r["n"]}
    # (c) contrôlé, quadratique
    Xq = X.copy(); Xq["dist2"] = d ** 2 / 100.0
    r = h.mco(Xq.values, y)
    out["controle_quadratique"] = {"b1": r["beta"][-2], "b2_100": r["beta"][-1], "se1": r["se_hc1"][-2], "se2": r["se_hc1"][-1], "r2": r["r2"],
                                   "gamma_0": -r["beta"][-2], "gamma_10": -(r["beta"][-2] + 2 * r["beta"][-1] * 10 / 100.0)}
    # (d) contrôlé, aire urbaine contiguë (d <= 16 km : Hull, Gatineau, Aylmer)
    m = d <= 16.0
    r = h.mco(X.values[m], y[m])
    out["controle_urbain_16km"] = {"gamma": -r["beta"][-1], "se": r["se_hc1"][-1], "r2": r["r2"], "n": int(m.sum())}
    # (e) modèle de référence (secteurs + max(d-3, 0))
    modele = modele or cas1.hedonique(v)
    c = modele["coef"]
    out["reference"] = {"gamma_au_dela_3km": -float(c.loc["dist_centre", "coef"]), "se": float(c.loc["dist_centre", "se_hc1"]),
                        "t": float(c.loc["dist_centre", "t"]), "r2": modele["r2"],
                        "secteurs": {s: float(c.loc[f"sect_{s}", "coef"]) for s in cas1.SECTEURS[1:]},
                        "secteurs_pct": {s: float(c.loc[f"sect_{s}", "effet_pct"]) for s in cas1.SECTEURS[1:]}}
    # (f) gradient implicite entre secteurs : médianes par secteur vs distance médiane
    gs = v.groupby("secteur").agg(dist=("dist_centre_km", "median"), prix=("prix_vente", "median"), n=("prix_vente", "size"))
    out["secteurs_medianes"] = gs
    return out

def gradient_exact_vs_exponentiel(P: dict = PARAMS) -> pd.DataFrame:
    d = np.array([0.0, 5.0, 10.0, 15.0, 20.0, 25.0])
    pe, pa = prix_logement(d, P), prix_exponentiel(d, P)
    return pd.DataFrame({"d_km": d, "p_exact": pe, "p_exponentiel": pa, "ecart_pct": pa / pe - 1, "q_m2": consommation_logement(d, P),
                         "rente_fonciere": rente_fonciere(d, P), "densite_relative": densite(d, P)})


# ----------------------------------------------------------------------------- cas 4 : décomposition interrégionale
def decomposition_rmr(pf: pd.DataFrame, reference: str = "Montréal") -> pd.DataFrame:
    """ln(V/n) = ln(RNE/n) - ln(R) par RMR, et écarts à la RMR de référence."""
    s = cas4.sommaire(pf)["par_rmr"].copy()
    s["rne_par_logement"] = s["rne"] / s["logements"]
    s["ln_V_n"] = np.log(s["valeur_par_logement"]); s["ln_RNE_n"] = np.log(s["rne_par_logement"]); s["ln_R"] = np.log(s["taux_pondere"])
    ref = s.loc[reference]
    s["ecart_ln_V_n"] = s["ln_V_n"] - ref["ln_V_n"]
    s["contrib_rne"] = s["ln_RNE_n"] - ref["ln_RNE_n"]
    s["contrib_taux"] = -(s["ln_R"] - ref["ln_R"])
    s["ratio_valeur"] = s["valeur_par_logement"] / ref["valeur_par_logement"]
    return s.sort_values("valeur_par_logement", ascending=False)


# ----------------------------------------------------------------------------- exemple intégral (impression)
def fmt(x, dec=0):
    s = f"{x:,.{dec}f}".replace(",", " ").replace(".", ",")
    return s

def main():
    P = PARAMS
    print("=== 5.A  Modèle monocentrique simulé (paramètres pédagogiques déclarés) ===")
    print({k: v for k, v in P.items()})
    g = gradient_exponentiel(P); dbar = frontiere_urbaine(P)
    print(f"gradient exponentiel au centre gamma = t/(beta y) = {g:.5f} par km ({g*100:.2f} %/km) ; demi-distance ln2/gamma = {demi_distance(g):.1f} km")
    print(f"frontière urbaine d_bar (r = ra) = {dbar:.2f} km ; p(d_bar) = {prix_logement(dbar, P):.2f} $/m²/an ; élasticité rente/prix 1/alpha = {elasticite_rente_prix(P['alpha']):.3f}")
    print(f"valeur du terrain au centre r0/i = {fmt(valeur_terrain(P['r0'], P['i']))} $/m² ; à la frontière ra/i = {valeur_terrain(P['ra'], P['i']):.0f} $/m²")
    tab = gradient_exact_vs_exponentiel(P)
    print(tab.round(3).to_string(index=False))
    print(f"q(0) = {consommation_logement(0.0, P):.1f} m² ; q(10) = {consommation_logement(10.0, P):.1f} m² ; q(d_bar) = {consommation_logement(dbar, P):.1f} m²")
    print("\n--- Alonso : frontières entre usages (km) ---")
    (R1, t1), (R2, t2), (R3, t3) = ALONSO.values()
    f12, f23 = frontiere_usages(R1, t1, R2, t2), frontiere_usages(R2, t2, R3, t3)
    f3a = (R3 - P["ra"]) / t3
    print(f"commerces/dense : {f12:.2f} km ; dense/pavillonnaire : {f23:.2f} km ; pavillonnaire/agricole : {f3a:.2f} km")
    print("\n--- von Thünen (exercice 5.1) ---")
    for k, q in THUNEN.items():
        print(f"{k}: R0 = {q['Y']*(q['p']-q['c']):.0f} $/ha/an ; pente = {q['Y']*q['f']:.0f} $/ha/an/km ; distance limite = {distance_limite_thunen(**q):.2f} km")
    a, b = THUNEN.values()
    print(f"frontière maraîchage/grandes cultures : {frontiere_usages(a['Y']*(a['p']-a['c']), a['Y']*a['f'], b['Y']*(b['p']-b['c']), b['Y']*b['f']):.2f} km")

    print("\n=== 5.B  Cas 1 (Gatineau, DONNÉES SYNTHÉTIQUES) : gradient prix/distance ===")
    d1 = cas1.charger(); v = d1["ventes"]; m = cas1.hedonique(v)
    tr = tranches_distance(v)
    print(tr.round(2).to_string())
    gr = gradients_cas1(v, m)
    for k in ["brut_prix_m2", "controle_lineaire", "controle_urbain_16km"]:
        r = gr[k]; print(f"{k}: gamma = {r['gamma']*100:.3f} %/km (e.-t. {r['se']*100:.3f}) ; R² = {r['r2']:.3f} ; n = {r['n']} ; demi-distance = {demi_distance(r['gamma']):.1f} km")
    q = gr["controle_quadratique"]
    print(f"controle_quadratique: b1 = {q['b1']*100:.3f} %/km (e.-t. {q['se1']*100:.3f}), b2/100 = {q['b2_100']*100:.4f} (e.-t. {q['se2']*100:.4f}) ; gamma(0) = {q['gamma_0']*100:.3f}, gamma(10 km) = {q['gamma_10']*100:.3f} %/km ; R² = {q['r2']:.3f}")
    rf = gr["reference"]
    print(f"reference (traite.cas1.hedonique) : dist_centre = -{rf['gamma_au_dela_3km']*100:.3f} %/km au-delà de 3 km (e.-t. {rf['se']*100:.3f}, t = {rf['t']:.2f}) ; R² = {rf['r2']:.3f}")
    print("  secteurs (réf. Hull), coef. et effet % :", {k: (round(v_, 4), round(rf['secteurs_pct'][k] * 100, 1)) for k, v_ in rf["secteurs"].items()})
    print(gr["secteurs_medianes"].round(2).to_string())
    gs = gr["secteurs_medianes"]
    pente = (math.log(gs.loc["Gatineau", "prix"]) - math.log(gs.loc["Hull", "prix"])) / (gs.loc["Gatineau", "dist"] - gs.loc["Hull", "dist"])
    print(f"gradient implicite Hull -> Gatineau (médianes) : {pente*100:.3f} %/km ; Hull -> Buckingham : {(math.log(gs.loc['Buckingham','prix'])-math.log(gs.loc['Hull','prix']))/(gs.loc['Buckingham','dist']-gs.loc['Hull','dist'])*100:.3f} %/km")
    print(f"sujet du cas 1 : dist_centre = {d1['sujet']['dist_centre_km']} km ; effet distance du modèle de référence = {max(d1['sujet']['dist_centre_km']-3,0)*rf['gamma_au_dela_3km']*100:.2f} %")
    # ce que vaut 1 km de plus pour le sujet selon la spécification contrôlée
    V_ind = 551029.0  # indication brute du ch. 16
    print(f"à partir de l'indication brute du ch. 16 ({fmt(V_ind)} $), 1 km de plus vaut {fmt(-V_ind*gr['controle_lineaire']['gamma'])} $ (spécification contrôlée) ; 5 km : {fmt(V_ind*(math.exp(-5*gr['controle_lineaire']['gamma'])-1))} $")

    print("\n=== 5.C  Cas 4 (portefeuille, DONNÉES SYNTHÉTIQUES) : valeur par logement par RMR ===")
    d4 = cas4.charger(); dec = decomposition_rmr(d4["portefeuille"])
    cols = ["n", "logements", "valeur_par_logement", "rne_par_logement", "taux_pondere", "ratio_valeur", "ecart_ln_V_n", "contrib_rne", "contrib_taux"]
    print(dec[cols].round(4).to_string())
    for rmr in dec.index:
        if rmr != "Montréal":
            e = dec.loc[rmr]; part = e["contrib_rne"] / e["ecart_ln_V_n"] if e["ecart_ln_V_n"] != 0 else float("nan")
            print(f"  {rmr}: écart ln = {e['ecart_ln_V_n']:.3f} dont RNE/log. {e['contrib_rne']:.3f} ({part*100:.0f} %) et taux {e['contrib_taux']:.3f} ({(1-part)*100:.0f} %)")

    print("\n=== 5.D  Rente résiduelle et levier ===")
    CN = 397567.33; dep = 0.30 * CN; V = 551029.0
    VL = valeur_residuelle(V, CN, dep)
    print(f"C_N (ch. 21) = {fmt(CN,2)} $ ; dépréciation hypothétique 30 % = {fmt(dep,2)} $ ; bâtiment déprécié = {fmt(CN-dep,2)} $")
    print(f"V (indication brute ch. 16) = {fmt(V)} $ -> V_L = {fmt(VL)} $ ({VL/V*100:.1f} % de V) ; levier V/V_L = {levier_residuel(V, VL):.2f}")
    for dv in (-0.05, 0.05, -0.10):
        VL2 = valeur_residuelle(V * (1 + dv), CN, dep)
        print(f"  V {dv:+.0%} -> V_L = {fmt(VL2)} $ ({VL2/VL-1:+.1%})")
    print(f"élasticité de la rente foncière au prix du logement (alpha = {P['alpha']}) : {elasticite_rente_prix(P['alpha']):.2f} ; +10 % de prix -> rente x {1.1**elasticite_rente_prix(P['alpha']):.3f}")
    cz = decomposition_capozza(P["ra"], P["i"], 25.0, rente_fonciere(10.0, P), 60.0)
    print("Capozza-Helsley à 10 km ($/m²) :", {k: round(v_, 1) for k, v_ in cz.items()})

    print("\n=== 5.E  Capitalisation de la taxe foncière ===")
    for k, t_ in TAUX_2026.items():
        print(f"  {k}: {t_:.5f} $/100 $ = {t_/100*10000:.1f} pb de la valeur")
    V0 = 500000.0; i = 0.05
    for pb in (10, 22.4, 50):
        c = capitalisation_taxe(V0, 0.0070, 0.0070 + pb / 10000, i)
        print(f"  écart de {pb} pb sur {fmt(V0)} $ : dT = {fmt(c['dT'])} $/an ; dV approché = {fmt(c['dV_approx'])} $ ; dV exact = {fmt(c['dV_exact'])} $ ({c['dV_pct_exact']*100:.2f} %) ; à 60 % : {fmt(0.6*c['dV_approx'])} $")
    ecart_qc_gat = (TAUX_2026["Québec (1 à 5 logements)"] - TAUX_2026["Gatineau (taux de base)"]) / 100
    print(f"écart nominal Québec (1-5 log.) - Gatineau (base) = {ecart_qc_gat*10000:.1f} pb ; sur {fmt(V0)} $ : {fmt(ecart_qc_gat*V0)} $/an ; capitalisé à 5 % : {fmt(ecart_qc_gat*V0/i)} $ ; à 60 % : {fmt(0.6*ecart_qc_gat*V0/i)} $")
    for ville, F in FACTEUR_2026.items():
        key = "Gatineau (taux de base)" if ville == "Gatineau" else "Montréal (taxe foncière générale, hors arrondissement)"
        print(f"  taux effectif {ville} : {TAUX_2026[key]:.4f}/{F} = {taux_effectif(TAUX_2026[key], F):.4f} $/100 $ de valeur marchande")
    print(f"  Québec avec un facteur hypothétique de 1,10 : {taux_effectif(TAUX_2026['Québec (1 à 5 logements)'], 1.10):.4f} $/100 $")
    sujet_role = d1["sujet"]["valeur_role_2025"]
    T_gat = sujet_role * TAUX_2026["Gatineau (taux de base)"] / 100; T_sco = sujet_role * TAUX_2026["Taxe scolaire 2026-2027 (toute la province)"] / 100
    print(f"sujet du cas 1 (rôle fictif {fmt(sujet_role)} $) : taxe générale au taux de base Gatineau = {fmt(T_gat)} $ ; scolaire (sans l'exemption de 25 000 $) = {fmt(T_sco)} $")
    print(f"  écart 6+ vs 1-5 logements à Québec : {(TAUX_2026['Québec (6 logements et plus)']-TAUX_2026['Québec (1 à 5 logements)'])*100:.1f} pb ; sur un immeuble de 2 000 000 $ : {fmt((TAUX_2026['Québec (6 logements et plus)']-TAUX_2026['Québec (1 à 5 logements)'])/100*2e6)} $/an ; capitalisé à 5,5 % : {fmt((TAUX_2026['Québec (6 logements et plus)']-TAUX_2026['Québec (1 à 5 logements)'])/100*2e6/0.055)} $")
    G = grille_capitalisation(V0, [5, 10, 25, 50, 100], [0.04, 0.05, 0.06, 0.07])
    print("grille (kappa = 1, V0 = 500 000 $) :\n", G.round(0).to_string())
    # coût d'usage : part de la taxe
    tau_eff = taux_effectif(TAUX_2026["Gatineau (taux de base)"], 1.02) / 100 + TAUX_2026["Taxe scolaire 2026-2027 (toute la province)"] / 100
    print(f"taux effectif total Gatineau (base + scolaire) ≈ {tau_eff*100:.4f} $/100 $ ; part de la taxe dans un coût d'usage de (i=4,12 %, tau, delta=1,5 %, g=2 %) : {tau_eff/(0.0412+tau_eff+0.015-0.02)*100:.1f} %")

if __name__ == "__main__":
    main()