#!/usr/bin/env python3
"""Chapitre 3 — Offre, demande et cycles : calculs de l'exemple intégral et des tableaux.

Trois blocs :
1. Modèle des quatre quadrants (DiPasquale–Wheaton) : équilibre stationnaire en forme fermée, statique comparative
   (taux de capitalisation, demande, coût de construction). PARAMÈTRES PÉDAGOGIQUES (marché type, données synthétiques).
2. Dynamique stock-flux trimestrielle avec délai de livraison et inoccupation naturelle : simulation d'un cycle
   après un choc de demande, racines caractéristiques du modèle linéaire (condition d'oscillation), sensibilité.
3. Données RÉELLES lues dans data/raw : mises en chantier au Québec (Statistique Canada 34-10-0135-01 et 34-10-0143-01),
   taux de la Banque du Canada (Valet : V39079, BD.CDN.5YR.DQ.YLD, V80691335, groupe A4), IPC Québec (18-10-0004-01),
   IPLN Québec (18-10-0205-01) ; coût d'usage du logement (Poterba) et ratio prix/loyer implicite ; mois d'inventaire
   (chiffres APCIQ T3 2026 cités dans sources/notes-E.md).

Exécuter depuis la racine : python3 analysis/python/cases/ch03_cycles.py  (imprime tous les chiffres cités au chapitre 3).
Testé par analysis/python/tests/test_ch03.py.
"""
from __future__ import annotations
import csv, math, pathlib, sys
from dataclasses import dataclass, replace
import numpy as np
import pandas as pd

ROOT = pathlib.Path(__file__).resolve().parents[3]
sys.path.insert(0, str(ROOT / "analysis/python"))
from traite import finance as fin  # noqa: E402  (P = R / i : capitalisation directe de la bibliothèque)

RAW = ROOT / "data/raw"
GRAINE = 20261031

# =============================================================================================================
# 1. Quatre quadrants — paramètres pédagogiques d'un marché type (SYNTHÉTIQUE)
# =============================================================================================================

@dataclass(frozen=True)
class Parametres4Q:
    """Marché type : S en logements, R en $/logement/an, P en $/logement, C en logements/an.
    Calibration pédagogique (ordres de grandeur) : loyer 1 200 $/mois ≈ loyer moyen SCHL RMR de Québec oct. 2025
    (1 227 $, tous logements) ; taux de capitalisation 5 % ≈ fourchette CBRE Québec T3 2026 (4,25-6,00 %) ;
    élasticité-prix de la demande d'espace −0,5 et taux de retrait δ = 1 %/an : hypothèses déclarées (non estimées)."""
    S0: float = 100_000.0          # stock initial (logements)
    R0: float = 14_400.0           # loyer annuel d'équilibre initial ($/logement/an) = 1 200 $/mois
    eta: float = -0.5              # élasticité-prix de la demande d'espace au point initial
    i: float = 0.05                # taux de capitalisation (marché des actifs)
    delta: float = 0.01            # taux de retrait du stock (démolitions, conversions), par an
    c0: float = 240_000.0          # seuil de coût : prix en deçà duquel rien ne se construit ($/logement)

    @property
    def b(self) -> float:
        """Pente de la demande D = a − bR ; b = −η S0 / R0."""
        return -self.eta * self.S0 / self.R0

    @property
    def a(self) -> float:
        return self.S0 + self.b * self.R0

    @property
    def C0(self) -> float:
        return self.delta * self.S0

    @property
    def P0(self) -> float:
        return fin.capitalisation_directe(self.R0, self.i)

    @property
    def c1(self) -> float:
        """Pente de l'offre de construction C = (P − c0)/c1, telle que C0 = δ S0 au prix initial."""
        return (self.P0 - self.c0) / self.C0


def equilibre_4q(p: Parametres4Q, a: float | None = None, i: float | None = None, c0: float | None = None) -> dict:
    """Équilibre stationnaire des quatre quadrants en forme fermée :
    S* = (a − b i c0) / (1 + b i c1 δ) ; R* = (a − S*)/b ; P* = R*/i ; C* = δ S*."""
    a = p.a if a is None else a
    i = p.i if i is None else i
    c0 = p.c0 if c0 is None else c0
    b, c1, d = p.b, p.c1, p.delta
    S = (a - b * i * c0) / (1.0 + b * i * c1 * d)
    R = (a - S) / b
    P = fin.capitalisation_directe(R, i)
    C = d * S
    assert abs(C - (P - c0) / c1) < 1e-6 * max(1.0, C)
    return {"S": S, "R": R, "P": P, "C": C, "a": a, "i": i, "c0": c0}


def equilibre_4q_iteratif(p: Parametres4Q, a: float | None = None, i: float | None = None, c0: float | None = None,
                          n_iter: int = 2000) -> dict:
    """Même équilibre par itération de point fixe S ← C(P(R(S)))/δ (vérification numérique de la forme fermée)."""
    a = p.a if a is None else a
    i = p.i if i is None else i
    c0 = p.c0 if c0 is None else c0
    S = p.S0
    for _ in range(n_iter):
        R = (a - S) / p.b
        P = fin.capitalisation_directe(R, i)
        C = max(0.0, (P - c0) / p.c1)
        S = 0.5 * S + 0.5 * C / p.delta       # relaxation pour garantir la convergence
    R = (a - S) / p.b
    return {"S": S, "R": R, "P": fin.capitalisation_directe(R, i), "C": p.delta * S}


def court_terme_4q(p: Parametres4Q, a: float | None = None, i: float | None = None) -> dict:
    """Court terme : stock fixe à S0 ; le loyer et le prix absorbent le choc."""
    a = p.a if a is None else a
    i = p.i if i is None else i
    R = (a - p.S0) / p.b
    P = fin.capitalisation_directe(R, i)
    return {"S": p.S0, "R": R, "P": P, "C": max(0.0, (P - p.c0) / p.c1)}


def statique_comparative(p: Parametres4Q) -> pd.DataFrame:
    """Trois chocs : baisse du taux de capitalisation à 4 %, hausse de la demande de 10 %, hausse du coût de 10 %."""
    base = equilibre_4q(p)
    chocs = {
        "Référence": dict(),
        "Taux i : 5 % → 4 %": dict(i=0.04),
        "Demande a : +10 %": dict(a=1.10 * p.a),
        "Coût c0 : +10 %": dict(c0=1.10 * p.c0),
    }
    rows = []
    for nom, kw in chocs.items():
        lt = equilibre_4q(p, **kw)
        ct = court_terme_4q(p, a=kw.get("a"), i=kw.get("i")) if "c0" not in kw else {"R": base["R"], "P": base["P"]}
        rows.append({"choc": nom, "R_ct": ct["R"], "P_ct": ct["P"], "S_lt": lt["S"], "R_lt": lt["R"], "P_lt": lt["P"], "C_lt": lt["C"],
                     "dS_pct": 100 * (lt["S"] / base["S"] - 1), "dR_ct_pct": 100 * (ct["R"] / base["R"] - 1),
                     "dR_lt_pct": 100 * (lt["R"] / base["R"] - 1), "dP_ct_pct": 100 * (ct["P"] / base["P"] - 1),
                     "dP_lt_pct": 100 * (lt["P"] / base["P"] - 1)})
    return pd.DataFrame(rows).set_index("choc")


def elasticites_4q(p: Parametres4Q) -> dict:
    """Élasticités au point initial : construction/prix (court terme) et stock stationnaire/demande (long terme)."""
    eps_C = (1.0 / p.c1) * (p.P0 / p.C0)                      # (dC/dP)(P/C)
    kappa = 1.0 / (p.b * p.i * p.c1)                          # réactivité de l'offre à l'écart de stock (par an)
    base = equilibre_4q(p)
    plus = equilibre_4q(p, a=1.01 * p.a)
    eps_S_a = (plus["S"] / base["S"] - 1) / 0.01               # élasticité du stock stationnaire à l'ordonnée de demande
    return {"eps_C_P": eps_C, "kappa": kappa, "eps_S_a": eps_S_a, "eta": p.eta}


# =============================================================================================================
# 2. Dynamique stock-flux (SYNTHÉTIQUE) — simulation trimestrielle et racines caractéristiques
# =============================================================================================================

@dataclass(frozen=True)
class ParametresCycle:
    p4q: Parametres4Q = Parametres4Q()
    v_nat: float = 0.03            # taux d'inoccupation naturel (SCHL RMR Québec oct. 2025 : 2,4 % ; Montréal 2,9 %)
    lam: float = 0.75              # vitesse d'ajustement du loyer : ΔR/R = −λ (v − v*) par trimestre
    L: int = 6                     # délai de livraison (trimestres) entre la mise en chantier et l'ajout au stock
    T: int = 100                   # horizon (trimestres)
    t_choc: int = 8                # trimestre où le choc commence
    n_choc: int = 12               # durée de montée du choc (trimestres) : la demande croît linéairement puis se stabilise
    choc_demande: float = 0.06     # hausse permanente de l'ordonnée de demande au terme de la montée (+6 %)
    mult_c1: float = 0.25          # multiplicateur de la pente d'offre c1 (0,25 : offre réactive, κ = 0,48 ; 1 : κ = 0,12)

    @property
    def c1(self) -> float:
        return self.p4q.c1 * self.mult_c1

    @property
    def c0(self) -> float:
        """Seuil de coût recalibré pour que l'état initial (S0, R0, P0, C0 = δ S0) soit stationnaire quelle que soit la pente c1."""
        return self.p4q.P0 - self.c1 * self.p4q.C0

    @property
    def kappa(self) -> float:
        """Réactivité annuelle de l'offre à l'écart de stock : κ = 1/(b i c1)."""
        return 1.0 / (self.p4q.b * self.p4q.i * self.c1)


def simuler_cycle(pc: ParametresCycle) -> pd.DataFrame:
    """Trimestriel. Demande occupée D = a − bR ; inoccupation v = 1 − D/S ; loyer R_{t+1} = R_t[1 − λ(v_t − v*)] ;
    prix P = R/i ; mises en chantier C = max(0, (P − c0)/c1)/4 ; stock S_{t+1} = (1 − δ/4) S_t + C_{t−L}.
    L'état initial est stationnaire avec v = v* (l'ordonnée a est ajustée pour que D0 = (1 − v*) S0)."""
    p = pc.p4q
    c1, c0 = pc.c1, pc.c0
    b = p.b
    a0 = (1.0 - pc.v_nat) * p.S0 + b * p.R0
    S, R = p.S0, p.R0
    starts = [p.C0 / 4.0] * (pc.L + 1)     # chantiers en cours au régime stationnaire
    rows = []
    for t in range(pc.T + 1):
        montee = min(1.0, max(0.0, (t - pc.t_choc + 1) / pc.n_choc)) if pc.n_choc > 0 else float(t >= pc.t_choc)
        a = a0 * (1.0 + pc.choc_demande * montee)
        D = min(S, a - b * R)
        v = 1.0 - D / S
        P = fin.capitalisation_directe(R, p.i)
        C = max(0.0, (P - c0) / c1) / 4.0
        rows.append({"t": t, "annee": t / 4.0, "S": S, "D": D, "v": v, "R": R, "P": P, "C": C, "a": a})
        starts.append(C)
        livraison = starts[-(pc.L + 1)]
        S = (1.0 - p.delta / 4.0) * S + livraison
        R = R * (1.0 - pc.lam * (v - pc.v_nat))
    df = pd.DataFrame(rows)
    df["R_mois"] = df["R"] / 12.0
    df["C_annuel"] = df["C"] * 4.0
    return df


def etat_stationnaire_cycle(pc: ParametresCycle, apres_choc: bool = True) -> dict:
    """État stationnaire du modèle dynamique (avec inoccupation naturelle) : v = v*, D = (1−v*)S, C = δS/4."""
    p = pc.p4q
    b = p.b
    a0 = (1.0 - pc.v_nat) * p.S0 + b * p.R0
    a = a0 * (1.0 + pc.choc_demande) if apres_choc else a0
    # (1−v*) S = a − bR ; R = i P ; P = c0 + c1 δ S  →  (1−v*) S = a − b i (c0 + c1 δ S)
    c1, c0 = pc.c1, pc.c0
    S = (a - b * p.i * c0) / ((1.0 - pc.v_nat) + b * p.i * c1 * p.delta)
    P = c0 + c1 * p.delta * S
    R = p.i * P
    return {"S": S, "R": R, "P": P, "C": p.delta * S, "v": pc.v_nat}


def resume_cycle(df: pd.DataFrame, pc: ParametresCycle) -> dict:
    """Mesures du cycle : pointe et creux de loyer (écart au stationnaire final), pointe d'inoccupation, période."""
    ss = etat_stationnaire_cycle(pc)
    apres = df[df["t"] >= pc.t_choc]
    R_rel = apres["R"] / ss["R"] - 1.0
    i_max = int(R_rel.idxmax())
    # creux après la première pointe
    post = R_rel.loc[i_max:]
    i_min = int(post.idxmin())
    # seconde pointe (pour la période) : maximum après le creux
    post2 = R_rel.loc[i_min:]
    i_max2 = int(post2.idxmax()) if len(post2) > 1 else i_max
    periode = (i_max2 - i_max) if i_max2 > i_max and post2.max() > 0 else float("nan")
    return {"R_pointe_pct": 100 * R_rel.max(), "t_pointe": df.loc[i_max, "t"] - pc.t_choc, "R_creux_pct": 100 * post.min(),
            "t_creux": df.loc[i_min, "t"] - pc.t_choc, "v_max_pct": 100 * apres["v"].max(), "v_min_pct": 100 * apres["v"].min(),
            "C_max_annuel": apres["C_annuel"].max(), "periode_trim": periode,
            "S_final_ecart_pct": 100 * (df["S"].iloc[-1] / ss["S"] - 1), "v_final_pct": 100 * df["v"].iloc[-1],
            "amortissement": (post2.max() / R_rel.max()) if periode == periode else float("nan")}


def sensibilite_cycle(pc: ParametresCycle, delais=(2, 4, 6, 8, 12), mult_c1=(0.25, 0.5, 1.0), lams=(0.5, 0.75, 1.5)) -> pd.DataFrame:
    rows = []
    for L in delais:
        for m in mult_c1:
            for lam in lams:
                r = resume_cycle(simuler_cycle(replace(pc, L=L, mult_c1=m, lam=lam)), replace(pc, L=L, mult_c1=m, lam=lam))
                rows.append({"L": L, "mult_c1": m, "lam": lam, "kappa": replace(pc, mult_c1=m).kappa, **r})
    return pd.DataFrame(rows)


def racines_caracteristiques(kappa: float, delta: float, L: int = 1) -> dict:
    """Modèle linéaire annuel S_{t+1} = (1−δ) S_t + κ (a − S_{t−L}) + cte, κ = 1/(b i c1).
    Polynôme : r^{L+1} − (1−δ) r^L + κ = 0. Retourne le module maximal, la nature des racines et la période."""
    coeffs = np.zeros(L + 2)
    coeffs[0] = 1.0
    coeffs[1] = -(1.0 - delta)
    coeffs[-1] = kappa
    r = np.roots(coeffs)
    dom = r[np.argmax(np.abs(r))]
    complexe = bool(abs(dom.imag) > 1e-9)
    periode = (2 * math.pi / abs(np.angle(dom))) if complexe else float("inf")
    return {"racines": r, "module": float(abs(dom)), "oscille": complexe, "periode": periode, "stable": bool(abs(dom) < 1.0)}


def condition_oscillation_L1(kappa: float, delta: float) -> dict:
    """Cas L = 1 en forme fermée : r² − (1−δ) r + κ = 0 ; racines complexes ssi κ > (1−δ)²/4 ; module √κ ;
    période 2π / arccos((1−δ)/(2√κ))."""
    seuil = (1.0 - delta) ** 2 / 4.0
    if kappa > seuil:
        return {"seuil": seuil, "oscille": True, "module": math.sqrt(kappa), "periode": 2 * math.pi / math.acos((1 - delta) / (2 * math.sqrt(kappa))),
                "stable": kappa < 1.0}
    disc = math.sqrt((1 - delta) ** 2 - 4 * kappa)
    r1, r2 = ((1 - delta) + disc) / 2, ((1 - delta) - disc) / 2
    return {"seuil": seuil, "oscille": False, "module": max(abs(r1), abs(r2)), "periode": float("inf"), "stable": max(abs(r1), abs(r2)) < 1.0}


# =============================================================================================================
# 3. Données RÉELLES
# =============================================================================================================

def _lire_statcan(fichier: pathlib.Path, filtre) -> list[dict]:
    with open(fichier, encoding="utf-8-sig") as f:
        return [r for r in csv.DictReader(f, delimiter=";") if filtre(r) and r["VALEUR"]]


def mises_en_chantier_quebec() -> pd.DataFrame:
    """34-10-0135-01 : Québec, toutes régions, trimestriel, non désaisonnalisé, mises en chantier par type d'unité.
    Périodes étiquetées par le premier mois du trimestre (1948-01 = T1 1948)."""
    rows = _lire_statcan(RAW / "statcan_34100135/34100135.csv",
                         lambda r: r["GÉO"] == "Québec" and r["Désaisonnalisation"] == "Non désaisonnalisés"
                         and r["Estimations de logement"] == "Logements mis en chantier")
    df = pd.DataFrame([{"periode": r["PÉRIODE DE RÉFÉRENCE"], "type": r["Type d'unité"], "unites": float(r["VALEUR"]), "vecteur": r["VECTEUR"]} for r in rows])
    df["annee"] = df["periode"].str[:4].astype(int)
    df["trimestre"] = df["periode"].str[5:7].map({"01": 1, "04": 2, "07": 3, "10": 4})
    df["date"] = pd.PeriodIndex(df["annee"].astype(str) + "Q" + df["trimestre"].astype(str), freq="Q").to_timestamp()
    return df.sort_values(["type", "date"]).reset_index(drop=True)


def mises_en_chantier_annuelles(df: pd.DataFrame | None = None) -> pd.DataFrame:
    """Sommes annuelles (années complètes à 4 trimestres) par type ; part des maisons individuelles et des collectifs."""
    df = mises_en_chantier_quebec() if df is None else df
    piv = df.pivot_table(index="annee", columns="type", values="unites", aggfunc="sum")
    n_trim = df[df["type"] == "Total d'unités"].groupby("annee")["trimestre"].count()
    piv["trimestres"] = n_trim
    piv = piv.rename(columns={"Total d'unités": "total", "Maisons individuelles": "individuelles", "Logements collectifs": "collectifs",
                              "Appartements et autres logements": "appartements", "Logements jumelés": "jumeles", "Logements en bande": "bande"})
    piv["part_individuelles"] = piv["individuelles"] / piv["total"]
    piv["part_collectifs"] = piv["collectifs"] / piv["total"]
    piv["part_appartements"] = piv["appartements"] / piv["total"]
    return piv


def mises_en_chantier_montreal_annuelles(geo: str = "Montréal, Québec") -> pd.DataFrame:
    """34-10-0143-01 : RMR de Montréal (centres de 10 000 hab. et plus), mensuel, total d'unités, sommes annuelles.
    La colonne « mois » compte les mois disponibles : dans le fichier du 2026-10-31, juillet à septembre 2025 sont ABSENTS
    pour Montréal, de sorte que la somme 2025 (21 138) ne couvre que 9 mois et n'est pas un total annuel."""
    rows = _lire_statcan(RAW / "statcan_34100143/34100143.csv",
                         lambda r: r["GÉO"] == geo and r["Estimations de logements"] == "Logements mis en chantier"
                         and r["Type d'unité"] == "Total d'unités")
    s = pd.Series({r["PÉRIODE DE RÉFÉRENCE"]: float(r["VALEUR"]) for r in rows})
    s.index = pd.to_datetime(s.index)
    ann = s.groupby(s.index.year).sum()
    mois = s.groupby(s.index.year).count()
    return pd.DataFrame({"total": ann, "mois": mois})


def _valet(nom: str) -> pd.Series:
    rows = list(csv.reader(open(RAW / "bdc_valet" / f"{nom}.csv", encoding="utf-8-sig")))
    i = next(k for k, r in enumerate(rows) if r and r[0] == "date")
    s = pd.Series({r[0]: float(r[1]) for r in rows[i + 1:] if len(r) > 1 and r[1]})
    s.index = pd.to_datetime(s.index)
    return s.sort_index()


def taux_bdc_mensuels() -> pd.DataFrame:
    """Moyennes mensuelles : taux hypothécaire affiché 5 ans (V80691335), obligation GdC 5 ans et 10 ans, taux directeur (V39079, depuis 2009 dans Valet)."""
    out = pd.DataFrame({
        "hyp5": _valet("V80691335").resample("MS").mean(),
        "gdc5": _valet("BD.CDN.5YR.DQ.YLD").resample("MS").mean(),
        "gdc10": _valet("BD.CDN.10YR.DQ.YLD").resample("MS").mean(),
        "directeur": _valet("V39079").resample("MS").mean(),
    })
    return out


def changements_taux_directeur(depuis: str = "2020-01-01") -> pd.DataFrame:
    """Dates d'effet (Valet, lendemain de l'annonce) et niveaux du taux cible du financement à un jour (V39079)."""
    s = _valet("V39079")
    ch = s[s.diff().fillna(0) != 0]
    ch = ch[ch.index >= depuis]
    prev = s.shift(1).reindex(ch.index)
    return pd.DataFrame({"avant": prev.values, "apres": ch.values, "variation_pb": ((ch - prev) * 100).round().astype(int).values}, index=ch.index)


def taux_effectifs_a4(mois: str = "2026-08-01") -> dict:
    """Groupe A4 (fonds avancés) : V122667783 = assuré, fixe 5 ans et plus ; V122667789 = non assuré, fixe 5 ans et plus
    (correspondance d'après les métadonnées du fichier ; valeurs d'août 2026 citées dans notes-C/E : 4,12 % et 4,58 %)."""
    rows = list(csv.reader(open(RAW / "bdc_valet/A4_RATES_MORTGAGES.csv", encoding="utf-8-sig")))
    i = next(k for k, r in enumerate(rows) if r and r[0] == "date")
    hdr = rows[i]
    rec = next(r for r in rows[i + 1:] if r and r[0] == mois)
    d = dict(zip(hdr, rec))
    # libellés des séries (bloc SERIES du fichier)
    labels = {r[0]: r[1] for r in rows[:i] if len(r) >= 2 and r[0].startswith("V1226677")}
    desc = {r[0]: (r[1], r[2]) for r in rows[:i] if len(r) >= 3 and r[0].startswith("V1226677")}
    def _cle(assure: bool):
        cible = "Residential mortgages, insured" if assure else "Residential mortgages, uninsured"
        return [k for k, (lab, de) in desc.items() if lab == "Fixed rate, 5 years and over" and de.startswith("Mortgages - Funds advanced") and cible in de]
    ka, kn = _cle(True), _cle(False)
    return {"mois": mois, "assure_fixe_5ans_plus": float(d[ka[0]]), "non_assure_fixe_5ans_plus": float(d[kn[0]]),
            "serie_assure": ka[0], "serie_non_assure": kn[0]}


def ipc_quebec(produit: str = "Ensemble") -> pd.Series:
    """18-10-0004-01, Québec, non désaisonnalisé, 2002 = 100 (Ensemble : v41691783 ; Loyer : v41691818)."""
    rows = _lire_statcan(RAW / "statcan_18100004/18100004.csv", lambda r: r["GÉO"] == "Québec" and r["Produits et groupes de produits"] == produit)
    s = pd.Series({r["PÉRIODE DE RÉFÉRENCE"] + "-01": float(r["VALEUR"]) for r in rows})
    s.index = pd.to_datetime(s.index)
    return s.sort_index()


def ipln_quebec(geo: str = "Québec", composante: str = "Total (maison et terrain)") -> pd.Series:
    rows = _lire_statcan(RAW / "statcan_18100205/18100205.csv", lambda r: r["GÉO"] == geo and r["Indices des prix des logements neufs"] == composante)
    s = pd.Series({r["PÉRIODE DE RÉFÉRENCE"] + "-01": float(r["VALEUR"]) for r in rows})
    s.index = pd.to_datetime(s.index)
    return s.sort_index()


# --- coût d'usage -------------------------------------------------------------------------------------------

TAUX_TAXE_QUEBEC_2026 = 0.007464 + 0.0007899   # Ville de Québec 1-5 log. (0,7464 $/100 $) + taxe scolaire 2026-2027 (0,07899 $/100 $)
DELTA_ENTRETIEN_RISQUE = 0.025                 # dépréciation + entretien + prime de risque : HYPOTHÈSE PÉDAGOGIQUE déclarée


def cout_usage(P: float, i: float, tau: float, delta: float, g: float) -> float:
    """u = P (i + τ + δ − g) : coût annuel d'occupation d'un logement de prix P (Poterba 1984)."""
    return P * (i + tau + delta - g)


def taux_cout_usage(i: float, tau: float, delta: float, g: float) -> float:
    return i + tau + delta - g


def ratio_prix_loyer(i: float, tau: float, delta: float, g: float) -> float:
    """À l'équilibre R = u  →  P/R = 1/(i + τ + δ − g) (multiple du loyer annuel)."""
    u = taux_cout_usage(i, tau, delta, g)
    if u <= 0:
        return float("inf")
    return 1.0 / u


def serie_cout_usage(tau: float = TAUX_TAXE_QUEBEC_2026, delta: float = DELTA_ENTRETIEN_RISQUE, debut: str = "2001-01-01") -> pd.DataFrame:
    """Série mensuelle réelle : i = taux hypothécaire affiché 5 ans ; π = inflation IPC Québec sur 12 mois ;
    g_ipc = π (anticipation = inflation courante) ; g_ipln = croissance annualisée de l'IPLN Québec sur 5 ans ;
    u et P/R pour chaque hypothèse. τ et δ constants (déclarés)."""
    bdc = taux_bdc_mensuels()
    ipc = ipc_quebec()
    pi = (ipc / ipc.shift(12) - 1.0) * 100
    ipln = ipln_quebec()
    g_ipln = ((ipln / ipln.shift(60)) ** (1 / 5) - 1.0) * 100
    df = pd.DataFrame({"hyp5": bdc["hyp5"], "gdc5": bdc["gdc5"], "directeur": bdc["directeur"], "pi": pi, "g_ipln": g_ipln}).loc[debut:]
    df = df.dropna(subset=["hyp5", "pi"])
    df["reel"] = df["hyp5"] - df["pi"]
    df["u_ipc"] = df["hyp5"] + 100 * (tau + delta) - df["pi"]
    df["u_ipln"] = df["hyp5"] + 100 * (tau + delta) - df["g_ipln"]
    df["pr_ipc"] = 100.0 / df["u_ipc"]
    df["pr_ipln"] = (100.0 / df["u_ipln"]).where(df["u_ipln"] > 0)
    return df


DATES_TABLEAU_COUT_USAGE = ["2001-01-01", "2008-01-01", "2012-01-01", "2020-07-01", "2023-10-01", "2026-09-01"]


def exemple_cout_usage(prix: float = 463_500.0, i_affiche: float = 0.0619, i_effectif: float = 0.0458, pi: float = 0.031,
                       g_ipln: float = 0.013, tau: float = TAUX_TAXE_QUEBEC_2026, delta: float = DELTA_ENTRETIEN_RISQUE) -> dict:
    """Exemple intégral : unifamiliale au prix médian APCIQ RMR de Québec T3 2026 (463 500 $) ; taux affiché 5 ans 6,19 %
    (BdC 2026-10-28) ; taux effectif non assuré fixe ≥ 5 ans 4,58 % (A4, août 2026) ; inflation IPC Québec 3,1 % ;
    IPLN Québec +1,3 % (sept. 2026). Semi-élasticité : d ln(P/R) / di = −1/u."""
    out = {}
    for nom, i in (("affiche", i_affiche), ("effectif", i_effectif)):
        for gnom, g in (("ipc", pi), ("ipln", g_ipln)):
            u = taux_cout_usage(i, tau, delta, g)
            out[f"u_{nom}_{gnom}"] = u
            out[f"cout_annuel_{nom}_{gnom}"] = cout_usage(prix, i, tau, delta, g)
            out[f"cout_mensuel_{nom}_{gnom}"] = cout_usage(prix, i, tau, delta, g) / 12
            out[f"pr_{nom}_{gnom}"] = ratio_prix_loyer(i, tau, delta, g)
    # sensibilité : +100 pb sur le taux effectif, anticipation = IPC
    u0 = taux_cout_usage(i_effectif, tau, delta, pi)
    u1 = taux_cout_usage(i_effectif + 0.01, tau, delta, pi)
    out["pr_effectif_ipc_plus100pb"] = 1 / u1
    out["var_pr_plus100pb_pct"] = 100 * ((1 / u1) / (1 / u0) - 1)
    out["semi_elasticite"] = -1 / u0
    # loyer de référence SCHL (RMR Québec oct. 2025, appartements 3 ch. et plus) pour comparaison
    out["loyer_schl_3ch_quebec"] = 1475.0
    out["loyer_schl_2ch_quebec"] = 1277.0
    out["prix_implicite_loyer_3ch_effectif_ipc"] = 1475.0 * 12 / u0
    return out


# --- indicateurs d'alerte (APCIQ T3 2026, chiffres cités dans sources/notes-E.md) --------------------------------

APCIQ_T3_2026 = {
    # territoire : (ventes T3 2026, inscriptions en vigueur, médiane unifamiliale)
    "Province": (20_935, 42_408, 500_000),
    "RMR Montréal": (9_265, 20_350, 645_500),
    "RMR Québec": (2_161, 2_204, 463_500),
    "RMR Gatineau": (1_044, 2_256, 480_000),
    "RMR Sherbrooke": (569, 825, 490_750),
    "RMR Trois-Rivières": (372, 447, 395_000),
    "RMR Saguenay": (371, 498, 365_500),
}


def mois_inventaire(inscriptions: float, ventes_trimestre: float) -> float:
    """Mois d'inventaire = inscriptions en vigueur / ventes mensuelles moyennes du trimestre."""
    return inscriptions / (ventes_trimestre / 3.0)


def indicateurs_apciq() -> pd.DataFrame:
    rows = []
    for terr, (v, insc, med) in APCIQ_T3_2026.items():
        rows.append({"territoire": terr, "ventes_t3": v, "inscriptions": insc, "ventes_mensuelles": v / 3.0,
                     "mois_inventaire": mois_inventaire(insc, v), "mediane_unifamiliale": med})
    return pd.DataFrame(rows).set_index("territoire")


# =============================================================================================================
# Exercices (chiffres des solutions, solutions/ch03.tex)
# =============================================================================================================

PARC_QUEBEC_2021 = 3_749_035   # ménages privés = logements privés occupés, Recensement 2021 (98-10-0240-01)


def exercices() -> dict:
    out = {}
    # ex. 1 : stock-flux avec δ = 0,5 %/an (hypothèse) et mises en chantier 2025
    C, S, d = 59_864.0, float(PARC_QUEBEC_2021), 0.005
    out["ex1_flux_pct"] = 100 * C / S
    out["ex1_retraits"] = d * S
    out["ex1_net"] = C - d * S
    out["ex1_net_pct"] = 100 * (C - d * S) / S
    out["ex1_annees_plus10"] = math.log(1.10) / math.log(1 + (C - d * S) / S)
    out["ex1_stock_stationnaire"] = C / d
    # ex. 2 : mois d'inventaire
    ind = indicateurs_apciq()
    out["ex2_quebec"] = ind.loc["RMR Québec", "mois_inventaire"]; out["ex2_montreal"] = ind.loc["RMR Montréal", "mois_inventaire"]
    out["ex2_gatineau"] = ind.loc["RMR Gatineau", "mois_inventaire"]
    out["ex2_montreal_ventes_moins10"] = mois_inventaire(20_350, 9_265 * 0.9)
    # ex. 3 : quatre quadrants (référence et c0 + 10 %)
    p = Parametres4Q(); out["ex3_base"] = equilibre_4q(p); out["ex3_c0_plus10"] = equilibre_4q(p, c0=1.10 * p.c0)
    # ex. 4 : baisse du taux à 4 %
    out["ex4_ct"] = court_terme_4q(p, i=0.04); out["ex4_lt"] = equilibre_4q(p, i=0.04)
    # ex. 5 : coût d'usage, unifamiliale RMR Montréal 645 500 $ (APCIQ T3 2026), i = 4,58 %, τ = 0,9 % (hypothèse), δ = 2,5 %
    P, i, tau, dlt = 645_500.0, 0.0458, 0.009, 0.025
    for g, nom in ((0.031, "ipc"), (0.013, "ipln")):
        out[f"ex5_u_{nom}"] = taux_cout_usage(i, tau, dlt, g); out[f"ex5_mensuel_{nom}"] = cout_usage(P, i, tau, dlt, g) / 12
        out[f"ex5_pr_{nom}"] = ratio_prix_loyer(i, tau, dlt, g)
    out["ex5_pr_plus100pb"] = ratio_prix_loyer(i + 0.01, tau, dlt, 0.031)
    out["ex5_var_pr_plus100pb_pct"] = 100 * (out["ex5_pr_plus100pb"] / out["ex5_pr_ipc"] - 1)
    out["ex5_loyer_schl_3ch_montreal"] = 1626.0
    out["ex5_prix_implicite_1626"] = 1626.0 * 12 / out["ex5_u_ipc"]
    # ex. 6 : variantes de la simulation
    pc = ParametresCycle()
    for L, m in ((12, 0.5), (2, 0.25), (6, 0.25)):
        r = resume_cycle(simuler_cycle(replace(pc, L=L, mult_c1=m)), replace(pc, L=L, mult_c1=m))
        out[f"ex6_L{L}_m{m}"] = {k: r[k] for k in ("R_pointe_pct", "R_creux_pct", "v_max_pct", "periode_trim", "amortissement")}
    # ex. 7 : condition d'oscillation
    out["ex7_k048"] = condition_oscillation_L1(0.48, p.delta); out["ex7_k030"] = condition_oscillation_L1(0.30, p.delta)
    out["ex7_c1_instable"] = 1.0 / (p.b * p.i)          # κ = 1 ⇔ c1 = 1/(b i)
    out["ex7_c1_seuil_oscillation"] = 1.0 / (p.b * p.i * (1 - p.delta) ** 2 / 4)
    return out


# =============================================================================================================
# Impression
# =============================================================================================================

def _f(x, d=0):
    return f"{x:,.{d}f}".replace(",", " ")


def main() -> None:
    p = Parametres4Q()
    print("=== 1. Quatre quadrants (marché type, SYNTHÉTIQUE) ===")
    print(f"b = {p.b:.4f} log./$ ; a = {_f(p.a)} ; P0 = {_f(p.P0)} $ ; C0 = {_f(p.C0)} ; c1 = {p.c1:.2f} $ par (log./an)")
    eq = equilibre_4q(p); it = equilibre_4q_iteratif(p)
    print(f"Équilibre fermé : S*={_f(eq['S'])} R*={_f(eq['R'])} P*={_f(eq['P'])} C*={_f(eq['C'])} | itératif S={_f(it['S'],1)}")
    el = elasticites_4q(p)
    print(f"Élasticités : construction/prix = {el['eps_C_P']:.2f} ; kappa = {el['kappa']:.4f} ; stock stationnaire/demande = {el['eps_S_a']:.3f}")
    sc = statique_comparative(p)
    pd.set_option("display.width", 200); pd.set_option("display.float_format", lambda x: f"{x:,.2f}")
    print(sc[["R_ct", "P_ct", "S_lt", "R_lt", "P_lt", "C_lt", "dS_pct", "dR_ct_pct", "dR_lt_pct", "dP_ct_pct", "dP_lt_pct"]])

    print("\n=== 2. Dynamique stock-flux (SYNTHÉTIQUE) ===")
    pc = ParametresCycle()
    df = simuler_cycle(pc); ss = etat_stationnaire_cycle(pc); r = resume_cycle(df, pc)
    print(f"kappa (référence dynamique) = {pc.kappa:.3f} ; c1 = {pc.c1:.1f} ; c0 recalibré = {_f(pc.c0)}")
    ss0 = etat_stationnaire_cycle(pc, apres_choc=False); print(f"Stationnaire initial : S={_f(ss0['S'])} R={_f(ss0['R'])} P={_f(ss0['P'])}")
    print(f"Stationnaire après choc : S={_f(ss['S'])} R={_f(ss['R'])} ({ss['R']/12:.0f} $/mois) P={_f(ss['P'])} C={_f(ss['C'])}/an")
    print({k: (round(v, 2) if isinstance(v, float) else v) for k, v in r.items()})
    for kap in (el["kappa"], pc.kappa):
        for L in (1, 2, 3):
            rc = racines_caracteristiques(kap, p.delta, L)
            print(f"Racines (kappa={kap:.3f}, L={L}) : module={rc['module']:.3f} oscille={rc['oscille']} période={rc['periode']:.1f} ans stable={rc['stable']}")
    for k in (0.12, 0.30, 0.60, 1.20):
        c = condition_oscillation_L1(k, p.delta)
        print(f"L=1, kappa={k}: seuil={c['seuil']:.4f} oscille={c['oscille']} module={c['module']:.3f} période={c['periode']:.2f} stable={c['stable']}")
    sens = sensibilite_cycle(pc)
    print(sens[(sens["lam"] == 0.75)][["L", "mult_c1", "kappa", "R_pointe_pct", "R_creux_pct", "v_max_pct", "C_max_annuel", "periode_trim", "amortissement"]].to_string(index=False))

    print("\n=== 3. Données réelles ===")
    ann = mises_en_chantier_annuelles()
    cols = ["total", "individuelles", "collectifs", "part_individuelles", "part_collectifs", "trimestres"]
    print(ann.loc[[1990, 2000, 2005, 2010, 2015, 2020, 2021, 2022, 2023, 2024, 2025, 2026], cols])
    mtl = mises_en_chantier_montreal_annuelles(); print("Montréal RMR (34-10-0143) :", mtl.loc[[2024, 2025, 2026]].to_dict())
    qc143 = mises_en_chantier_montreal_annuelles("Québec"); print("Québec centres 10 000+ (34-10-0143) :", qc143.loc[[2024, 2025, 2026]].to_dict())
    ch = changements_taux_directeur(); print(ch)
    print(taux_effectifs_a4())
    cu = serie_cout_usage()
    print(cu.loc[DATES_TABLEAU_COUT_USAGE, ["hyp5", "gdc5", "pi", "g_ipln", "reel", "u_ipc", "pr_ipc", "u_ipln", "pr_ipln"]])
    print("Dernier mois :", cu.index[-1].date(), cu.iloc[-1][["hyp5", "pi", "u_ipc", "pr_ipc"]].to_dict())
    print("IPC Québec sept. 2026 :", ipc_quebec().iloc[-1], "; a/a =", round((ipc_quebec().iloc[-1] / ipc_quebec().iloc[-13] - 1) * 100, 2))
    print("IPLN Québec sept. 2026 :", ipln_quebec().iloc[-1], "; a/a =", round((ipln_quebec().iloc[-1] / ipln_quebec().iloc[-13] - 1) * 100, 2))
    ex = exemple_cout_usage()
    for k, v in ex.items():
        print(f"  {k}: {v:,.4f}" if isinstance(v, float) else f"  {k}: {v}")
    print(indicateurs_apciq())
    print("\n=== Exercices ===")
    for k, v in exercices().items():
        print(f"  {k}: {v}")


if __name__ == "__main__":
    main()