spb/traite-evaluation-qc

1#!/usr/bin/env python32"""Chapitre 3 — Offre, demande et cycles : calculs de l'exemple intégral et des tableaux.34Trois blocs :51. Modèle des quatre quadrants (DiPasquale–Wheaton) : équilibre stationnaire en forme fermée, statique comparative6   (taux de capitalisation, demande, coût de construction). PARAMÈTRES PÉDAGOGIQUES (marché type, données synthétiques).72. Dynamique stock-flux trimestrielle avec délai de livraison et inoccupation naturelle : simulation d'un cycle8   après un choc de demande, racines caractéristiques du modèle linéaire (condition d'oscillation), sensibilité.93. 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),10   taux de la Banque du Canada (Valet : V39079, BD.CDN.5YR.DQ.YLD, V80691335, groupe A4), IPC Québec (18-10-0004-01),11   IPLN Québec (18-10-0205-01) ; coût d'usage du logement (Poterba) et ratio prix/loyer implicite ; mois d'inventaire12   (chiffres APCIQ T3 2026 cités dans sources/notes-E.md).1314Exécuter depuis la racine : python3 analysis/python/cases/ch03_cycles.py  (imprime tous les chiffres cités au chapitre 3).15Testé par analysis/python/tests/test_ch03.py.16"""17from __future__ import annotations18import csv, math, pathlib, sys19from dataclasses import dataclass, replace20import numpy as np21import pandas as pd2223ROOT = pathlib.Path(__file__).resolve().parents[3]24sys.path.insert(0, str(ROOT / "analysis/python"))25from traite import finance as fin  # noqa: E402  (P = R / i : capitalisation directe de la bibliothèque)2627RAW = ROOT / "data/raw"28GRAINE = 202610312930# =============================================================================================================31# 1. Quatre quadrants — paramètres pédagogiques d'un marché type (SYNTHÉTIQUE)32# =============================================================================================================3334@dataclass(frozen=True)35class Parametres4Q:36    """Marché type : S en logements, R en $/logement/an, P en $/logement, C en logements/an.37    Calibration pédagogique (ordres de grandeur) : loyer 1 200 $/mois ≈ loyer moyen SCHL RMR de Québec oct. 202538    (1 227 $, tous logements) ; taux de capitalisation 5 % ≈ fourchette CBRE Québec T3 2026 (4,25-6,00 %) ;39    élasticité-prix de la demande d'espace −0,5 et taux de retrait δ = 1 %/an : hypothèses déclarées (non estimées)."""40    S0: float = 100_000.0          # stock initial (logements)41    R0: float = 14_400.0           # loyer annuel d'équilibre initial ($/logement/an) = 1 200 $/mois42    eta: float = -0.5              # élasticité-prix de la demande d'espace au point initial43    i: float = 0.05                # taux de capitalisation (marché des actifs)44    delta: float = 0.01            # taux de retrait du stock (démolitions, conversions), par an45    c0: float = 240_000.0          # seuil de coût : prix en deçà duquel rien ne se construit ($/logement)4647    @property48    def b(self) -> float:49        """Pente de la demande D = a − bR ; b = −η S0 / R0."""50        return -self.eta * self.S0 / self.R05152    @property53    def a(self) -> float:54        return self.S0 + self.b * self.R05556    @property57    def C0(self) -> float:58        return self.delta * self.S05960    @property61    def P0(self) -> float:62        return fin.capitalisation_directe(self.R0, self.i)6364    @property65    def c1(self) -> float:66        """Pente de l'offre de construction C = (P − c0)/c1, telle que C0 = δ S0 au prix initial."""67        return (self.P0 - self.c0) / self.C0686970def equilibre_4q(p: Parametres4Q, a: float | None = None, i: float | None = None, c0: float | None = None) -> dict:71    """Équilibre stationnaire des quatre quadrants en forme fermée :72    S* = (a − b i c0) / (1 + b i c1 δ) ; R* = (a − S*)/b ; P* = R*/i ; C* = δ S*."""73    a = p.a if a is None else a74    i = p.i if i is None else i75    c0 = p.c0 if c0 is None else c076    b, c1, d = p.b, p.c1, p.delta77    S = (a - b * i * c0) / (1.0 + b * i * c1 * d)78    R = (a - S) / b79    P = fin.capitalisation_directe(R, i)80    C = d * S81    assert abs(C - (P - c0) / c1) < 1e-6 * max(1.0, C)82    return {"S": S, "R": R, "P": P, "C": C, "a": a, "i": i, "c0": c0}838485def equilibre_4q_iteratif(p: Parametres4Q, a: float | None = None, i: float | None = None, c0: float | None = None,86                          n_iter: int = 2000) -> dict:87    """Même équilibre par itération de point fixe S ← C(P(R(S)))/δ (vérification numérique de la forme fermée)."""88    a = p.a if a is None else a89    i = p.i if i is None else i90    c0 = p.c0 if c0 is None else c091    S = p.S092    for _ in range(n_iter):93        R = (a - S) / p.b94        P = fin.capitalisation_directe(R, i)95        C = max(0.0, (P - c0) / p.c1)96        S = 0.5 * S + 0.5 * C / p.delta       # relaxation pour garantir la convergence97    R = (a - S) / p.b98    return {"S": S, "R": R, "P": fin.capitalisation_directe(R, i), "C": p.delta * S}99100101def court_terme_4q(p: Parametres4Q, a: float | None = None, i: float | None = None) -> dict:102    """Court terme : stock fixe à S0 ; le loyer et le prix absorbent le choc."""103    a = p.a if a is None else a104    i = p.i if i is None else i105    R = (a - p.S0) / p.b106    P = fin.capitalisation_directe(R, i)107    return {"S": p.S0, "R": R, "P": P, "C": max(0.0, (P - p.c0) / p.c1)}108109110def statique_comparative(p: Parametres4Q) -> pd.DataFrame:111    """Trois chocs : baisse du taux de capitalisation à 4 %, hausse de la demande de 10 %, hausse du coût de 10 %."""112    base = equilibre_4q(p)113    chocs = {114        "Référence": dict(),115        "Taux i : 5 % → 4 %": dict(i=0.04),116        "Demande a : +10 %": dict(a=1.10 * p.a),117        "Coût c0 : +10 %": dict(c0=1.10 * p.c0),118    }119    rows = []120    for nom, kw in chocs.items():121        lt = equilibre_4q(p, **kw)122        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"]}123        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"],124                     "dS_pct": 100 * (lt["S"] / base["S"] - 1), "dR_ct_pct": 100 * (ct["R"] / base["R"] - 1),125                     "dR_lt_pct": 100 * (lt["R"] / base["R"] - 1), "dP_ct_pct": 100 * (ct["P"] / base["P"] - 1),126                     "dP_lt_pct": 100 * (lt["P"] / base["P"] - 1)})127    return pd.DataFrame(rows).set_index("choc")128129130def elasticites_4q(p: Parametres4Q) -> dict:131    """Élasticités au point initial : construction/prix (court terme) et stock stationnaire/demande (long terme)."""132    eps_C = (1.0 / p.c1) * (p.P0 / p.C0)                      # (dC/dP)(P/C)133    kappa = 1.0 / (p.b * p.i * p.c1)                          # réactivité de l'offre à l'écart de stock (par an)134    base = equilibre_4q(p)135    plus = equilibre_4q(p, a=1.01 * p.a)136    eps_S_a = (plus["S"] / base["S"] - 1) / 0.01               # élasticité du stock stationnaire à l'ordonnée de demande137    return {"eps_C_P": eps_C, "kappa": kappa, "eps_S_a": eps_S_a, "eta": p.eta}138139140# =============================================================================================================141# 2. Dynamique stock-flux (SYNTHÉTIQUE) — simulation trimestrielle et racines caractéristiques142# =============================================================================================================143144@dataclass(frozen=True)145class ParametresCycle:146    p4q: Parametres4Q = Parametres4Q()147    v_nat: float = 0.03            # taux d'inoccupation naturel (SCHL RMR Québec oct. 2025 : 2,4 % ; Montréal 2,9 %)148    lam: float = 0.75              # vitesse d'ajustement du loyer : ΔR/R = −λ (v − v*) par trimestre149    L: int = 6                     # délai de livraison (trimestres) entre la mise en chantier et l'ajout au stock150    T: int = 100                   # horizon (trimestres)151    t_choc: int = 8                # trimestre où le choc commence152    n_choc: int = 12               # durée de montée du choc (trimestres) : la demande croît linéairement puis se stabilise153    choc_demande: float = 0.06     # hausse permanente de l'ordonnée de demande au terme de la montée (+6 %)154    mult_c1: float = 0.25          # multiplicateur de la pente d'offre c1 (0,25 : offre réactive, κ = 0,48 ; 1 : κ = 0,12)155156    @property157    def c1(self) -> float:158        return self.p4q.c1 * self.mult_c1159160    @property161    def c0(self) -> float:162        """Seuil de coût recalibré pour que l'état initial (S0, R0, P0, C0 = δ S0) soit stationnaire quelle que soit la pente c1."""163        return self.p4q.P0 - self.c1 * self.p4q.C0164165    @property166    def kappa(self) -> float:167        """Réactivité annuelle de l'offre à l'écart de stock : κ = 1/(b i c1)."""168        return 1.0 / (self.p4q.b * self.p4q.i * self.c1)169170171def simuler_cycle(pc: ParametresCycle) -> pd.DataFrame:172    """Trimestriel. Demande occupée D = a − bR ; inoccupation v = 1 − D/S ; loyer R_{t+1} = R_t[1 − λ(v_t − v*)] ;173    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}.174    L'état initial est stationnaire avec v = v* (l'ordonnée a est ajustée pour que D0 = (1 − v*) S0)."""175    p = pc.p4q176    c1, c0 = pc.c1, pc.c0177    b = p.b178    a0 = (1.0 - pc.v_nat) * p.S0 + b * p.R0179    S, R = p.S0, p.R0180    starts = [p.C0 / 4.0] * (pc.L + 1)     # chantiers en cours au régime stationnaire181    rows = []182    for t in range(pc.T + 1):183        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)184        a = a0 * (1.0 + pc.choc_demande * montee)185        D = min(S, a - b * R)186        v = 1.0 - D / S187        P = fin.capitalisation_directe(R, p.i)188        C = max(0.0, (P - c0) / c1) / 4.0189        rows.append({"t": t, "annee": t / 4.0, "S": S, "D": D, "v": v, "R": R, "P": P, "C": C, "a": a})190        starts.append(C)191        livraison = starts[-(pc.L + 1)]192        S = (1.0 - p.delta / 4.0) * S + livraison193        R = R * (1.0 - pc.lam * (v - pc.v_nat))194    df = pd.DataFrame(rows)195    df["R_mois"] = df["R"] / 12.0196    df["C_annuel"] = df["C"] * 4.0197    return df198199200def etat_stationnaire_cycle(pc: ParametresCycle, apres_choc: bool = True) -> dict:201    """État stationnaire du modèle dynamique (avec inoccupation naturelle) : v = v*, D = (1−v*)S, C = δS/4."""202    p = pc.p4q203    b = p.b204    a0 = (1.0 - pc.v_nat) * p.S0 + b * p.R0205    a = a0 * (1.0 + pc.choc_demande) if apres_choc else a0206    # (1−v*) S = a − bR ; R = i P ; P = c0 + c1 δ S  →  (1−v*) S = a − b i (c0 + c1 δ S)207    c1, c0 = pc.c1, pc.c0208    S = (a - b * p.i * c0) / ((1.0 - pc.v_nat) + b * p.i * c1 * p.delta)209    P = c0 + c1 * p.delta * S210    R = p.i * P211    return {"S": S, "R": R, "P": P, "C": p.delta * S, "v": pc.v_nat}212213214def resume_cycle(df: pd.DataFrame, pc: ParametresCycle) -> dict:215    """Mesures du cycle : pointe et creux de loyer (écart au stationnaire final), pointe d'inoccupation, période."""216    ss = etat_stationnaire_cycle(pc)217    apres = df[df["t"] >= pc.t_choc]218    R_rel = apres["R"] / ss["R"] - 1.0219    i_max = int(R_rel.idxmax())220    # creux après la première pointe221    post = R_rel.loc[i_max:]222    i_min = int(post.idxmin())223    # seconde pointe (pour la période) : maximum après le creux224    post2 = R_rel.loc[i_min:]225    i_max2 = int(post2.idxmax()) if len(post2) > 1 else i_max226    periode = (i_max2 - i_max) if i_max2 > i_max and post2.max() > 0 else float("nan")227    return {"R_pointe_pct": 100 * R_rel.max(), "t_pointe": df.loc[i_max, "t"] - pc.t_choc, "R_creux_pct": 100 * post.min(),228            "t_creux": df.loc[i_min, "t"] - pc.t_choc, "v_max_pct": 100 * apres["v"].max(), "v_min_pct": 100 * apres["v"].min(),229            "C_max_annuel": apres["C_annuel"].max(), "periode_trim": periode,230            "S_final_ecart_pct": 100 * (df["S"].iloc[-1] / ss["S"] - 1), "v_final_pct": 100 * df["v"].iloc[-1],231            "amortissement": (post2.max() / R_rel.max()) if periode == periode else float("nan")}232233234def 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:235    rows = []236    for L in delais:237        for m in mult_c1:238            for lam in lams:239                r = resume_cycle(simuler_cycle(replace(pc, L=L, mult_c1=m, lam=lam)), replace(pc, L=L, mult_c1=m, lam=lam))240                rows.append({"L": L, "mult_c1": m, "lam": lam, "kappa": replace(pc, mult_c1=m).kappa, **r})241    return pd.DataFrame(rows)242243244def racines_caracteristiques(kappa: float, delta: float, L: int = 1) -> dict:245    """Modèle linéaire annuel S_{t+1} = (1−δ) S_t + κ (a − S_{t−L}) + cte, κ = 1/(b i c1).246    Polynôme : r^{L+1} − (1−δ) r^L + κ = 0. Retourne le module maximal, la nature des racines et la période."""247    coeffs = np.zeros(L + 2)248    coeffs[0] = 1.0249    coeffs[1] = -(1.0 - delta)250    coeffs[-1] = kappa251    r = np.roots(coeffs)252    dom = r[np.argmax(np.abs(r))]253    complexe = bool(abs(dom.imag) > 1e-9)254    periode = (2 * math.pi / abs(np.angle(dom))) if complexe else float("inf")255    return {"racines": r, "module": float(abs(dom)), "oscille": complexe, "periode": periode, "stable": bool(abs(dom) < 1.0)}256257258def condition_oscillation_L1(kappa: float, delta: float) -> dict:259    """Cas L = 1 en forme fermée : r² − (1−δ) r + κ = 0 ; racines complexes ssi κ > (1−δ)²/4 ; module √κ ;260    période 2π / arccos((1−δ)/(2√κ))."""261    seuil = (1.0 - delta) ** 2 / 4.0262    if kappa > seuil:263        return {"seuil": seuil, "oscille": True, "module": math.sqrt(kappa), "periode": 2 * math.pi / math.acos((1 - delta) / (2 * math.sqrt(kappa))),264                "stable": kappa < 1.0}265    disc = math.sqrt((1 - delta) ** 2 - 4 * kappa)266    r1, r2 = ((1 - delta) + disc) / 2, ((1 - delta) - disc) / 2267    return {"seuil": seuil, "oscille": False, "module": max(abs(r1), abs(r2)), "periode": float("inf"), "stable": max(abs(r1), abs(r2)) < 1.0}268269270# =============================================================================================================271# 3. Données RÉELLES272# =============================================================================================================273274def _lire_statcan(fichier: pathlib.Path, filtre) -> list[dict]:275    with open(fichier, encoding="utf-8-sig") as f:276        return [r for r in csv.DictReader(f, delimiter=";") if filtre(r) and r["VALEUR"]]277278279def mises_en_chantier_quebec() -> pd.DataFrame:280    """34-10-0135-01 : Québec, toutes régions, trimestriel, non désaisonnalisé, mises en chantier par type d'unité.281    Périodes étiquetées par le premier mois du trimestre (1948-01 = T1 1948)."""282    rows = _lire_statcan(RAW / "statcan_34100135/34100135.csv",283                         lambda r: r["GÉO"] == "Québec" and r["Désaisonnalisation"] == "Non désaisonnalisés"284                         and r["Estimations de logement"] == "Logements mis en chantier")285    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])286    df["annee"] = df["periode"].str[:4].astype(int)287    df["trimestre"] = df["periode"].str[5:7].map({"01": 1, "04": 2, "07": 3, "10": 4})288    df["date"] = pd.PeriodIndex(df["annee"].astype(str) + "Q" + df["trimestre"].astype(str), freq="Q").to_timestamp()289    return df.sort_values(["type", "date"]).reset_index(drop=True)290291292def mises_en_chantier_annuelles(df: pd.DataFrame | None = None) -> pd.DataFrame:293    """Sommes annuelles (années complètes à 4 trimestres) par type ; part des maisons individuelles et des collectifs."""294    df = mises_en_chantier_quebec() if df is None else df295    piv = df.pivot_table(index="annee", columns="type", values="unites", aggfunc="sum")296    n_trim = df[df["type"] == "Total d'unités"].groupby("annee")["trimestre"].count()297    piv["trimestres"] = n_trim298    piv = piv.rename(columns={"Total d'unités": "total", "Maisons individuelles": "individuelles", "Logements collectifs": "collectifs",299                              "Appartements et autres logements": "appartements", "Logements jumelés": "jumeles", "Logements en bande": "bande"})300    piv["part_individuelles"] = piv["individuelles"] / piv["total"]301    piv["part_collectifs"] = piv["collectifs"] / piv["total"]302    piv["part_appartements"] = piv["appartements"] / piv["total"]303    return piv304305306def mises_en_chantier_montreal_annuelles(geo: str = "Montréal, Québec") -> pd.DataFrame:307    """34-10-0143-01 : RMR de Montréal (centres de 10 000 hab. et plus), mensuel, total d'unités, sommes annuelles.308    La colonne « mois » compte les mois disponibles : dans le fichier du 2026-10-31, juillet à septembre 2025 sont ABSENTS309    pour Montréal, de sorte que la somme 2025 (21 138) ne couvre que 9 mois et n'est pas un total annuel."""310    rows = _lire_statcan(RAW / "statcan_34100143/34100143.csv",311                         lambda r: r["GÉO"] == geo and r["Estimations de logements"] == "Logements mis en chantier"312                         and r["Type d'unité"] == "Total d'unités")313    s = pd.Series({r["PÉRIODE DE RÉFÉRENCE"]: float(r["VALEUR"]) for r in rows})314    s.index = pd.to_datetime(s.index)315    ann = s.groupby(s.index.year).sum()316    mois = s.groupby(s.index.year).count()317    return pd.DataFrame({"total": ann, "mois": mois})318319320def _valet(nom: str) -> pd.Series:321    rows = list(csv.reader(open(RAW / "bdc_valet" / f"{nom}.csv", encoding="utf-8-sig")))322    i = next(k for k, r in enumerate(rows) if r and r[0] == "date")323    s = pd.Series({r[0]: float(r[1]) for r in rows[i + 1:] if len(r) > 1 and r[1]})324    s.index = pd.to_datetime(s.index)325    return s.sort_index()326327328def taux_bdc_mensuels() -> pd.DataFrame:329    """Moyennes mensuelles : taux hypothécaire affiché 5 ans (V80691335), obligation GdC 5 ans et 10 ans, taux directeur (V39079, depuis 2009 dans Valet)."""330    out = pd.DataFrame({331        "hyp5": _valet("V80691335").resample("MS").mean(),332        "gdc5": _valet("BD.CDN.5YR.DQ.YLD").resample("MS").mean(),333        "gdc10": _valet("BD.CDN.10YR.DQ.YLD").resample("MS").mean(),334        "directeur": _valet("V39079").resample("MS").mean(),335    })336    return out337338339def changements_taux_directeur(depuis: str = "2020-01-01") -> pd.DataFrame:340    """Dates d'effet (Valet, lendemain de l'annonce) et niveaux du taux cible du financement à un jour (V39079)."""341    s = _valet("V39079")342    ch = s[s.diff().fillna(0) != 0]343    ch = ch[ch.index >= depuis]344    prev = s.shift(1).reindex(ch.index)345    return pd.DataFrame({"avant": prev.values, "apres": ch.values, "variation_pb": ((ch - prev) * 100).round().astype(int).values}, index=ch.index)346347348def taux_effectifs_a4(mois: str = "2026-08-01") -> dict:349    """Groupe A4 (fonds avancés) : V122667783 = assuré, fixe 5 ans et plus ; V122667789 = non assuré, fixe 5 ans et plus350    (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 %)."""351    rows = list(csv.reader(open(RAW / "bdc_valet/A4_RATES_MORTGAGES.csv", encoding="utf-8-sig")))352    i = next(k for k, r in enumerate(rows) if r and r[0] == "date")353    hdr = rows[i]354    rec = next(r for r in rows[i + 1:] if r and r[0] == mois)355    d = dict(zip(hdr, rec))356    # libellés des séries (bloc SERIES du fichier)357    labels = {r[0]: r[1] for r in rows[:i] if len(r) >= 2 and r[0].startswith("V1226677")}358    desc = {r[0]: (r[1], r[2]) for r in rows[:i] if len(r) >= 3 and r[0].startswith("V1226677")}359    def _cle(assure: bool):360        cible = "Residential mortgages, insured" if assure else "Residential mortgages, uninsured"361        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]362    ka, kn = _cle(True), _cle(False)363    return {"mois": mois, "assure_fixe_5ans_plus": float(d[ka[0]]), "non_assure_fixe_5ans_plus": float(d[kn[0]]),364            "serie_assure": ka[0], "serie_non_assure": kn[0]}365366367def ipc_quebec(produit: str = "Ensemble") -> pd.Series:368    """18-10-0004-01, Québec, non désaisonnalisé, 2002 = 100 (Ensemble : v41691783 ; Loyer : v41691818)."""369    rows = _lire_statcan(RAW / "statcan_18100004/18100004.csv", lambda r: r["GÉO"] == "Québec" and r["Produits et groupes de produits"] == produit)370    s = pd.Series({r["PÉRIODE DE RÉFÉRENCE"] + "-01": float(r["VALEUR"]) for r in rows})371    s.index = pd.to_datetime(s.index)372    return s.sort_index()373374375def ipln_quebec(geo: str = "Québec", composante: str = "Total (maison et terrain)") -> pd.Series:376    rows = _lire_statcan(RAW / "statcan_18100205/18100205.csv", lambda r: r["GÉO"] == geo and r["Indices des prix des logements neufs"] == composante)377    s = pd.Series({r["PÉRIODE DE RÉFÉRENCE"] + "-01": float(r["VALEUR"]) for r in rows})378    s.index = pd.to_datetime(s.index)379    return s.sort_index()380381382# --- coût d'usage -------------------------------------------------------------------------------------------383384TAUX_TAXE_QUEBEC_2026 = 0.007464 + 0.0007899   # Ville de Québec 1-5 log. (0,7464 $/100 $) + taxe scolaire 2026-2027 (0,07899 $/100 $)385DELTA_ENTRETIEN_RISQUE = 0.025                 # dépréciation + entretien + prime de risque : HYPOTHÈSE PÉDAGOGIQUE déclarée386387388def cout_usage(P: float, i: float, tau: float, delta: float, g: float) -> float:389    """u = P (i + τ + δ − g) : coût annuel d'occupation d'un logement de prix P (Poterba 1984)."""390    return P * (i + tau + delta - g)391392393def taux_cout_usage(i: float, tau: float, delta: float, g: float) -> float:394    return i + tau + delta - g395396397def ratio_prix_loyer(i: float, tau: float, delta: float, g: float) -> float:398    """À l'équilibre R = u  →  P/R = 1/(i + τ + δ − g) (multiple du loyer annuel)."""399    u = taux_cout_usage(i, tau, delta, g)400    if u <= 0:401        return float("inf")402    return 1.0 / u403404405def serie_cout_usage(tau: float = TAUX_TAXE_QUEBEC_2026, delta: float = DELTA_ENTRETIEN_RISQUE, debut: str = "2001-01-01") -> pd.DataFrame:406    """Série mensuelle réelle : i = taux hypothécaire affiché 5 ans ; π = inflation IPC Québec sur 12 mois ;407    g_ipc = π (anticipation = inflation courante) ; g_ipln = croissance annualisée de l'IPLN Québec sur 5 ans ;408    u et P/R pour chaque hypothèse. τ et δ constants (déclarés)."""409    bdc = taux_bdc_mensuels()410    ipc = ipc_quebec()411    pi = (ipc / ipc.shift(12) - 1.0) * 100412    ipln = ipln_quebec()413    g_ipln = ((ipln / ipln.shift(60)) ** (1 / 5) - 1.0) * 100414    df = pd.DataFrame({"hyp5": bdc["hyp5"], "gdc5": bdc["gdc5"], "directeur": bdc["directeur"], "pi": pi, "g_ipln": g_ipln}).loc[debut:]415    df = df.dropna(subset=["hyp5", "pi"])416    df["reel"] = df["hyp5"] - df["pi"]417    df["u_ipc"] = df["hyp5"] + 100 * (tau + delta) - df["pi"]418    df["u_ipln"] = df["hyp5"] + 100 * (tau + delta) - df["g_ipln"]419    df["pr_ipc"] = 100.0 / df["u_ipc"]420    df["pr_ipln"] = (100.0 / df["u_ipln"]).where(df["u_ipln"] > 0)421    return df422423424DATES_TABLEAU_COUT_USAGE = ["2001-01-01", "2008-01-01", "2012-01-01", "2020-07-01", "2023-10-01", "2026-09-01"]425426427def exemple_cout_usage(prix: float = 463_500.0, i_affiche: float = 0.0619, i_effectif: float = 0.0458, pi: float = 0.031,428                       g_ipln: float = 0.013, tau: float = TAUX_TAXE_QUEBEC_2026, delta: float = DELTA_ENTRETIEN_RISQUE) -> dict:429    """Exemple intégral : unifamiliale au prix médian APCIQ RMR de Québec T3 2026 (463 500 $) ; taux affiché 5 ans 6,19 %430    (BdC 2026-10-28) ; taux effectif non assuré fixe ≥ 5 ans 4,58 % (A4, août 2026) ; inflation IPC Québec 3,1 % ;431    IPLN Québec +1,3 % (sept. 2026). Semi-élasticité : d ln(P/R) / di = −1/u."""432    out = {}433    for nom, i in (("affiche", i_affiche), ("effectif", i_effectif)):434        for gnom, g in (("ipc", pi), ("ipln", g_ipln)):435            u = taux_cout_usage(i, tau, delta, g)436            out[f"u_{nom}_{gnom}"] = u437            out[f"cout_annuel_{nom}_{gnom}"] = cout_usage(prix, i, tau, delta, g)438            out[f"cout_mensuel_{nom}_{gnom}"] = cout_usage(prix, i, tau, delta, g) / 12439            out[f"pr_{nom}_{gnom}"] = ratio_prix_loyer(i, tau, delta, g)440    # sensibilité : +100 pb sur le taux effectif, anticipation = IPC441    u0 = taux_cout_usage(i_effectif, tau, delta, pi)442    u1 = taux_cout_usage(i_effectif + 0.01, tau, delta, pi)443    out["pr_effectif_ipc_plus100pb"] = 1 / u1444    out["var_pr_plus100pb_pct"] = 100 * ((1 / u1) / (1 / u0) - 1)445    out["semi_elasticite"] = -1 / u0446    # loyer de référence SCHL (RMR Québec oct. 2025, appartements 3 ch. et plus) pour comparaison447    out["loyer_schl_3ch_quebec"] = 1475.0448    out["loyer_schl_2ch_quebec"] = 1277.0449    out["prix_implicite_loyer_3ch_effectif_ipc"] = 1475.0 * 12 / u0450    return out451452453# --- indicateurs d'alerte (APCIQ T3 2026, chiffres cités dans sources/notes-E.md) --------------------------------454455APCIQ_T3_2026 = {456    # territoire : (ventes T3 2026, inscriptions en vigueur, médiane unifamiliale)457    "Province": (20_935, 42_408, 500_000),458    "RMR Montréal": (9_265, 20_350, 645_500),459    "RMR Québec": (2_161, 2_204, 463_500),460    "RMR Gatineau": (1_044, 2_256, 480_000),461    "RMR Sherbrooke": (569, 825, 490_750),462    "RMR Trois-Rivières": (372, 447, 395_000),463    "RMR Saguenay": (371, 498, 365_500),464}465466467def mois_inventaire(inscriptions: float, ventes_trimestre: float) -> float:468    """Mois d'inventaire = inscriptions en vigueur / ventes mensuelles moyennes du trimestre."""469    return inscriptions / (ventes_trimestre / 3.0)470471472def indicateurs_apciq() -> pd.DataFrame:473    rows = []474    for terr, (v, insc, med) in APCIQ_T3_2026.items():475        rows.append({"territoire": terr, "ventes_t3": v, "inscriptions": insc, "ventes_mensuelles": v / 3.0,476                     "mois_inventaire": mois_inventaire(insc, v), "mediane_unifamiliale": med})477    return pd.DataFrame(rows).set_index("territoire")478479480# =============================================================================================================481# Exercices (chiffres des solutions, solutions/ch03.tex)482# =============================================================================================================483484PARC_QUEBEC_2021 = 3_749_035   # ménages privés = logements privés occupés, Recensement 2021 (98-10-0240-01)485486487def exercices() -> dict:488    out = {}489    # ex. 1 : stock-flux avec δ = 0,5 %/an (hypothèse) et mises en chantier 2025490    C, S, d = 59_864.0, float(PARC_QUEBEC_2021), 0.005491    out["ex1_flux_pct"] = 100 * C / S492    out["ex1_retraits"] = d * S493    out["ex1_net"] = C - d * S494    out["ex1_net_pct"] = 100 * (C - d * S) / S495    out["ex1_annees_plus10"] = math.log(1.10) / math.log(1 + (C - d * S) / S)496    out["ex1_stock_stationnaire"] = C / d497    # ex. 2 : mois d'inventaire498    ind = indicateurs_apciq()499    out["ex2_quebec"] = ind.loc["RMR Québec", "mois_inventaire"]; out["ex2_montreal"] = ind.loc["RMR Montréal", "mois_inventaire"]500    out["ex2_gatineau"] = ind.loc["RMR Gatineau", "mois_inventaire"]501    out["ex2_montreal_ventes_moins10"] = mois_inventaire(20_350, 9_265 * 0.9)502    # ex. 3 : quatre quadrants (référence et c0 + 10 %)503    p = Parametres4Q(); out["ex3_base"] = equilibre_4q(p); out["ex3_c0_plus10"] = equilibre_4q(p, c0=1.10 * p.c0)504    # ex. 4 : baisse du taux à 4 %505    out["ex4_ct"] = court_terme_4q(p, i=0.04); out["ex4_lt"] = equilibre_4q(p, i=0.04)506    # ex. 5 : coût d'usage, unifamiliale RMR Montréal 645 500 $ (APCIQ T3 2026), i = 4,58 %, τ = 0,9 % (hypothèse), δ = 2,5 %507    P, i, tau, dlt = 645_500.0, 0.0458, 0.009, 0.025508    for g, nom in ((0.031, "ipc"), (0.013, "ipln")):509        out[f"ex5_u_{nom}"] = taux_cout_usage(i, tau, dlt, g); out[f"ex5_mensuel_{nom}"] = cout_usage(P, i, tau, dlt, g) / 12510        out[f"ex5_pr_{nom}"] = ratio_prix_loyer(i, tau, dlt, g)511    out["ex5_pr_plus100pb"] = ratio_prix_loyer(i + 0.01, tau, dlt, 0.031)512    out["ex5_var_pr_plus100pb_pct"] = 100 * (out["ex5_pr_plus100pb"] / out["ex5_pr_ipc"] - 1)513    out["ex5_loyer_schl_3ch_montreal"] = 1626.0514    out["ex5_prix_implicite_1626"] = 1626.0 * 12 / out["ex5_u_ipc"]515    # ex. 6 : variantes de la simulation516    pc = ParametresCycle()517    for L, m in ((12, 0.5), (2, 0.25), (6, 0.25)):518        r = resume_cycle(simuler_cycle(replace(pc, L=L, mult_c1=m)), replace(pc, L=L, mult_c1=m))519        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")}520    # ex. 7 : condition d'oscillation521    out["ex7_k048"] = condition_oscillation_L1(0.48, p.delta); out["ex7_k030"] = condition_oscillation_L1(0.30, p.delta)522    out["ex7_c1_instable"] = 1.0 / (p.b * p.i)          # κ = 1 ⇔ c1 = 1/(b i)523    out["ex7_c1_seuil_oscillation"] = 1.0 / (p.b * p.i * (1 - p.delta) ** 2 / 4)524    return out525526527# =============================================================================================================528# Impression529# =============================================================================================================530531def _f(x, d=0):532    return f"{x:,.{d}f}".replace(",", " ")533534535def main() -> None:536    p = Parametres4Q()537    print("=== 1. Quatre quadrants (marché type, SYNTHÉTIQUE) ===")538    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)")539    eq = equilibre_4q(p); it = equilibre_4q_iteratif(p)540    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)}")541    el = elasticites_4q(p)542    print(f"Élasticités : construction/prix = {el['eps_C_P']:.2f} ; kappa = {el['kappa']:.4f} ; stock stationnaire/demande = {el['eps_S_a']:.3f}")543    sc = statique_comparative(p)544    pd.set_option("display.width", 200); pd.set_option("display.float_format", lambda x: f"{x:,.2f}")545    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"]])546547    print("\n=== 2. Dynamique stock-flux (SYNTHÉTIQUE) ===")548    pc = ParametresCycle()549    df = simuler_cycle(pc); ss = etat_stationnaire_cycle(pc); r = resume_cycle(df, pc)550    print(f"kappa (référence dynamique) = {pc.kappa:.3f} ; c1 = {pc.c1:.1f} ; c0 recalibré = {_f(pc.c0)}")551    ss0 = etat_stationnaire_cycle(pc, apres_choc=False); print(f"Stationnaire initial : S={_f(ss0['S'])} R={_f(ss0['R'])} P={_f(ss0['P'])}")552    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")553    print({k: (round(v, 2) if isinstance(v, float) else v) for k, v in r.items()})554    for kap in (el["kappa"], pc.kappa):555        for L in (1, 2, 3):556            rc = racines_caracteristiques(kap, p.delta, L)557            print(f"Racines (kappa={kap:.3f}, L={L}) : module={rc['module']:.3f} oscille={rc['oscille']} période={rc['periode']:.1f} ans stable={rc['stable']}")558    for k in (0.12, 0.30, 0.60, 1.20):559        c = condition_oscillation_L1(k, p.delta)560        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']}")561    sens = sensibilite_cycle(pc)562    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))563564    print("\n=== 3. Données réelles ===")565    ann = mises_en_chantier_annuelles()566    cols = ["total", "individuelles", "collectifs", "part_individuelles", "part_collectifs", "trimestres"]567    print(ann.loc[[1990, 2000, 2005, 2010, 2015, 2020, 2021, 2022, 2023, 2024, 2025, 2026], cols])568    mtl = mises_en_chantier_montreal_annuelles(); print("Montréal RMR (34-10-0143) :", mtl.loc[[2024, 2025, 2026]].to_dict())569    qc143 = mises_en_chantier_montreal_annuelles("Québec"); print("Québec centres 10 000+ (34-10-0143) :", qc143.loc[[2024, 2025, 2026]].to_dict())570    ch = changements_taux_directeur(); print(ch)571    print(taux_effectifs_a4())572    cu = serie_cout_usage()573    print(cu.loc[DATES_TABLEAU_COUT_USAGE, ["hyp5", "gdc5", "pi", "g_ipln", "reel", "u_ipc", "pr_ipc", "u_ipln", "pr_ipln"]])574    print("Dernier mois :", cu.index[-1].date(), cu.iloc[-1][["hyp5", "pi", "u_ipc", "pr_ipc"]].to_dict())575    print("IPC Québec sept. 2026 :", ipc_quebec().iloc[-1], "; a/a =", round((ipc_quebec().iloc[-1] / ipc_quebec().iloc[-13] - 1) * 100, 2))576    print("IPLN Québec sept. 2026 :", ipln_quebec().iloc[-1], "; a/a =", round((ipln_quebec().iloc[-1] / ipln_quebec().iloc[-13] - 1) * 100, 2))577    ex = exemple_cout_usage()578    for k, v in ex.items():579        print(f"  {k}: {v:,.4f}" if isinstance(v, float) else f"  {k}: {v}")580    print(indicateurs_apciq())581    print("\n=== Exercices ===")582    for k, v in exercices().items():583        print(f"  {k}: {v}")584585586if __name__ == "__main__":587    main()588