Pipeline de validation empirique de la méthode des comparables (Vrai Prix)
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
13 changed files +1,422 −0
added
.gitignore
+5 −0
@@ -0,0 +1,5 @@ | ||
| 1 | +out/ | |
| 2 | +data/ | |
| 3 | +*.pyc | |
| 4 | +__pycache__/ | |
| 5 | +.DS_Store | |
added
README.md
+33 −0
@@ -0,0 +1,33 @@ | ||
| 1 | +# Validation empirique de la méthode des comparables (Vrai Prix) | |
| 2 | + | |
| 3 | +Projet de recherche : rétro-évaluation à grande échelle de la méthode de comparaison des ventes sur les transactions | |
| 4 | +résidentielles unifamiliales du Québec (Vrai Prix, 2021-01 → 2026-09). Protocole : `docs/PROTOCOLE.md`. | |
| 5 | +Rapport de résultats : `docs/RAPPORT.md` / `docs/RAPPORT.pdf`. | |
| 6 | + | |
| 7 | +## Organisation | |
| 8 | + | |
| 9 | +``` | |
| 10 | +pipeline/ 01_dataset → 02_index → 03_models → 04_backtest → 05_analyze → 06_report | |
| 11 | +docs/ PROTOCOLE.md, RAPPORT.md(.pdf), header.tex | |
| 12 | +out/<version>/ sales.parquet, index_asof.npy, models_*.json, est_E*.npy, targets.parquet, tables/, fig/, results.json | |
| 13 | +data/ (nœud seulement) vraiprix-20261010.db — instantané de la base de production Vrai Prix | |
| 14 | +ref/ MUN.csv (répertoire des municipalités MAMH : code → MRC → région) | |
| 15 | +``` | |
| 16 | + | |
| 17 | +## Exécution | |
| 18 | + | |
| 19 | +Le calcul se fait sur **M2U64** (`~/comparables-validation`, venv `~/uqo-eval-terrains/env`, Python 3.14 : pandas, numpy, scipy, | |
| 20 | +lightgbm, matplotlib). Le laptop garde le code, les docs et une copie des sorties. | |
| 21 | + | |
| 22 | +```bash | |
| 23 | +rsync -a --exclude out --exclude data --exclude .git ~/Desktop/Validation-Comparables/ M2U64:comparables-validation/ | |
| 24 | +ssh M2U64 'cd ~/comparables-validation && CV_VERSION=2026-10-10 ./run_all.sh' # ~ 30 min | |
| 25 | +./make_pdf.sh # rapatrie out/ + docs/ et compile le PDF | |
| 26 | +``` | |
| 27 | + | |
| 28 | +Variables : `CV_VERSION` (dossier de sortie), `CV_SAMPLE=3000` (échantillon de cibles pour un essai), `CV_PROCS` (processus, 20 par défaut). | |
| 29 | + | |
| 30 | +## Données | |
| 31 | + | |
| 32 | +Instantané `vraiprix.db` pris sur M3U96a (`~/apps/vrai-prix/data`) le 2026-10-10 : 766 596 transactions, 3,75 M d'unités du rôle. | |
| 33 | +Jeu propre : ventes `unifamilial` jumelées au rôle, règles R0-R11 journalisées dans `out/<v>/exclusions.csv`. | |
added
docs/PROTOCOLE.md
+101 −0
@@ -0,0 +1,101 @@ | ||
| 1 | +# Revisiter la méthode de comparaison des ventes : protocole de validation empirique à grande échelle | |
| 2 | + | |
| 3 | +**Titre de travail (article A)** — *Revisiting the Sales Comparison Approach: A Large-Scale Empirical Validation of Comparable Selection, Price Adjustments, and Value Reconciliation* | |
| 4 | + | |
| 5 | +**Variante** — *Does the Sales Comparison Approach Actually Work? Evidence from Large-Scale Residential Property Transactions* | |
| 6 | + | |
| 7 | +Données : transactions historiques **Vrai Prix** (registre foncier du Québec jumelé au rôle d'évaluation), instantané du 2026-10-10. | |
| 8 | + | |
| 9 | +## 1. Question de recherche | |
| 10 | + | |
| 11 | +La méthode des comparables, telle qu'enseignée et pratiquée par les évaluateurs agréés, produit-elle réellement les estimations les plus fiables ? Et surtout, quelles règles de **sélection**, d'**ajustement**, de **pondération** et de **réconciliation** des comparables sont empiriquement justifiées ? | |
| 12 | + | |
| 13 | +| Dimension | Question testée | Où dans le pipeline | | |
| 14 | +|---|---|---| | |
| 15 | +| Sélection | Quels comparables sont réellement les plus pertinents (proximité, similarité physique, « moins ajustés ») ? | E1 (facteur `sel`), E2 | | |
| 16 | +| Nombre | Faut-il 1, 3, 5, 10, 20 ou 50 ventes ? | E1 (facteur `k`) | | |
| 17 | +| Ajustements | Quels ajustements sont statistiquement justifiés (grille conventionnelle, hédonique, locale, non linéaire, ratio au rôle) ? | E2 (facteur `adj`) | | |
| 18 | +| Pondération | Égalitaire, distance, récence, similarité ? | E2 (facteur `w`) | | |
| 19 | +| Réconciliation | Médiane, moyenne, moyenne pondérée, combinaison optimisée (M9) ? | E2 (facteur `recon`), M9 | | |
| 20 | +| Fiabilité | Dans quels marchés (densité, région, segment de prix, hétérogénéité) la méthode tient-elle ou défaille-t-elle ? | strates, H2-H5 | | |
| 21 | + | |
| 22 | +## 2. Protocole expérimental (rétro-évaluation « point dans le temps ») | |
| 23 | + | |
| 24 | +Pour chaque transaction cible : | |
| 25 | + | |
| 26 | +1. **Masquage** du prix de vente réel. | |
| 27 | +2. **Reconstitution du marché à la date d'évaluation** : seules les ventes conclues *strictement avant* la date de la cible sont admissibles (le bien lui-même, s'il a déjà été vendu, est exclu du bassin : on évalue par comparables, pas par vente répétée). L'indice temporel, les coefficients hédoniques et les modèles d'apprentissage sont ré-estimés avec les seules ventes antérieures à la période d'évaluation (indice : ré-estimé chaque mois avec toutes les ventes des mois précédents ; modèles : estimés sur les ventes < 2024-01 pour la validation 2024 et < 2025-01 pour le test 2025-26). | |
| 28 | +3. **Application** de chaque configuration (sélection × rayon × fenêtre × nombre × ajustement × pondération × réconciliation) sur exactement les mêmes cibles. | |
| 29 | +4. **Comparaison** avec le prix réel : erreur absolue en pourcentage, biais, part des estimations à ±5/10/20 %, couverture (part des cibles pour lesquelles la configuration produit une valeur), dispersion des prix ajustés (intervalle de valeur). | |
| 30 | + | |
| 31 | +### Séparation chronologique | |
| 32 | + | |
| 33 | +| Période | Rôle | Ventes | | |
| 34 | +|---|---|---| | |
| 35 | +| 2021-01 → 2023-12 | développement (estimation des grilles d'ajustement et des modèles) | ~ 200 000 | | |
| 36 | +| 2024 | validation : choix des configurations et des hyperparamètres | ~ 65 000 cibles | | |
| 37 | +| 2025-01 → 2026-09 | évaluation finale hors échantillon | ~ 125 000 cibles | | |
| 38 | + | |
| 39 | +Les comparables d'une cible de 2025 peuvent être des ventes de 2024 ou de 2025 antérieures à sa date : c'est l'information dont disposait un évaluateur ce jour-là (nous ignorons le délai de publication au registre foncier, ce qui avantage légèrement toutes les méthodes de la même façon). | |
| 40 | + | |
| 41 | +## 3. Expériences | |
| 42 | + | |
| 43 | +### E1 — surface empirique de performance Erreur(k, d, t) | |
| 44 | + | |
| 45 | +- k ∈ {1, 3, 5, 7, 10, 15, 20, 30, 50} comparables ; | |
| 46 | +- d ∈ {0,5 ; 1 ; 2 ; 5 ; 10 ; 25} km (rayon maximal) ; | |
| 47 | +- t ∈ {3, 6, 12, 24} mois (ancienneté maximale) ; | |
| 48 | +- sélection ∈ {`geo` les plus proches, `sim` les plus semblables physiquement, `mix` compromis} ; | |
| 49 | +- ajustement ∈ {`none` temps seulement, `hedo` grille hédonique} ; réconciliation = médiane. | |
| 50 | + | |
| 51 | +1 296 configurations × toutes les cibles. | |
| 52 | + | |
| 53 | +### E2 — plan factoriel des règles (trois ancres (k, d, t) = (3, 2 km, 6 mois), (5, 5 km, 12 mois), (10, 10 km, 24 mois)) | |
| 54 | + | |
| 55 | +- sélection ∈ {`geo`, `sim`, `mix`, `adj` = comparables nécessitant le moins d'ajustements} ; | |
| 56 | +- ajustement ∈ {`raw` aucun, `none` temps seul, `conv` grille conventionnelle stylisée, `hedo` grille hédonique globale (within-municipalité), `hedo_reg` grille régionale rétrécie, `local` régression sur le voisinage (≤ 10 km, ≥ 40 ventes, ridge vers la grille régionale), `gbm` ajustements non linéaires (LightGBM), `role` ratio prix/rôle (comparables de la même municipalité et du même rôle)} ; | |
| 57 | +- pondération ∈ {égale, 1/(0,25 + d), exp(−mois/6), 1/(0,1 + ajustement brut), produit des trois} ; | |
| 58 | +- réconciliation ∈ {médiane (pondérée), moyenne (pondérée) en dollars}. | |
| 59 | + | |
| 60 | +960 configurations. | |
| 61 | + | |
| 62 | +### Compétition de modèles | |
| 63 | + | |
| 64 | +| Modèle | Définition opérationnelle | | |
| 65 | +|---|---| | |
| 66 | +| M1 | médiane des 3 ventes les plus proches (≤ 2 km, ≤ 6 mois), temps seul | | |
| 67 | +| M2 | moyenne des 5 plus proches (≤ 5 km, ≤ 12 mois), temps seul | | |
| 68 | +| M3 | 5 comparables choisis par similarité physique | | |
| 69 | +| M4 | 5 plus proches pondérés par la distance | | |
| 70 | +| M5 | 5 plus proches pondérés par la récence | | |
| 71 | +| M6 | grille d'ajustement conventionnelle stylisée (½ du prix au m² pour l'aire, ½ de la valeur foncière au m² pour le terrain plafonnée à 15 %, −0,5 %/an d'âge plafonné à 20 %, ±5 % détaché/attaché) | | |
| 72 | +| M7 | grille d'ajustement hédonique (coefficients within-municipalité estimés sur les ventes antérieures) | | |
| 73 | +| M8a / M8b | grille locale (régression sur le voisinage) / ajustements non linéaires (GBM) | | |
| 74 | +| M9 | réconciliation optimisée : combinaison convexe (NNLS) de sept estimations par comparables, poids appris en 2024 | | |
| 75 | +| M10 | modèle hédonique global direct (β within + effets région, MRC, municipalité, unité de voisinage) | | |
| 76 | +| M11a / M11b | gradient boosting sans / avec la valeur au rôle | | |
| 77 | +| M12 | hybride : combinaison convexe de M11b et des estimations par comparables ajustées (GBM, hédonique) | | |
| 78 | +| M13 | ratio au rôle : rôle de la cible × ratio prix/rôle médian des 5 plus proches du même rôle | | |
| 79 | +| M0 | meilleure configuration de E1 ∪ E2 choisie sur la validation 2024, évaluée en 2025-26 | | |
| 80 | + | |
| 81 | +### Hypothèses | |
| 82 | + | |
| 83 | +- **H1** — un comparable très proche mais peu semblable peut être moins informatif qu'un comparable plus éloigné mais très semblable. | |
| 84 | +- **H2** — l'ajout de comparables améliore la précision jusqu'à un seuil, au-delà duquel l'hétérogénéité domine ; le seuil dépend de la densité. | |
| 85 | +- **H3** — les ajustements empiriques locaux surpassent les grilles fixes mais deviennent instables dans les petits échantillons. | |
| 86 | +- **H4** — la pondération optimale dépend de la liquidité du sous-marché. | |
| 87 | +- **H5** — les méthodes simples et transparentes rivalisent avec l'apprentissage automatique dans les marchés homogènes. | |
| 88 | + | |
| 89 | +### Inférence | |
| 90 | + | |
| 91 | +Différences d'erreurs appariées (mêmes cibles) avec intervalles de confiance par **bootstrap par grappes** (municipalité × trimestre), qui respecte les dépendances spatiales et temporelles. | |
| 92 | + | |
| 93 | +## 4. Ce que les données permettent — et ne permettent pas | |
| 94 | + | |
| 95 | +Caractéristiques observées pour chaque vente (registre foncier ⋈ rôle) : date, prix, coordonnées, municipalité, unité de voisinage, aire habitable (étages), superficie du terrain, année de construction, nombre d'étages, lien physique (détaché, jumelé, en rangée), genre de construction, valeur au rôle (en vigueur à la vente) et sa part terrain. **Non observés** : garage, piscine, nombre de chambres et de salles de bains, qualité et état, rénovations. Les ajustements correspondants ne peuvent donc pas être testés ici ; les erreurs résiduelles les contiennent. | |
| 96 | + | |
| 97 | +Nous validons la **précision prédictive** de règles professionnelles *simulées*, pas le jugement réel d'un évaluateur (il faudrait ses dossiers, ses comparables et ses ajustements). La distinction est maintenue partout dans le rapport. | |
| 98 | + | |
| 99 | +## 5. Article B (suite envisagée) | |
| 100 | + | |
| 101 | +Vers une méthode des comparables optimisée empiriquement : algorithme transparent de sélection/ajustement entraîné sur les ventes historiques et validé sur les ventes futures ; règles conditionnelles au type de bien, à la densité et à la qualité des données ; intégration dans UQO Éval pour montrer aux étudiants la comparaison entre règles traditionnelles et règles estimées. | |
added
docs/header.tex
+28 −0
@@ -0,0 +1,28 @@ | ||
| 1 | +% En-tête LaTeX — Grille empirique d'ajustement des terrains (UQO Éval) | |
| 2 | +\usepackage{fontspec} | |
| 3 | +\setmainfont{Helvetica Neue}[Ligatures=TeX] | |
| 4 | +\setsansfont{Helvetica Neue} | |
| 5 | +\setmonofont{Menlo}[Scale=0.82] | |
| 6 | +\usepackage{newunicodechar} | |
| 7 | +\newunicodechar{→}{\ensuremath{\rightarrow}}\newunicodechar{≥}{\ensuremath{\geq}}\newunicodechar{≤}{\ensuremath{\leq}}\newunicodechar{×}{\ensuremath{\times}}\newunicodechar{−}{\ensuremath{-}}\newunicodechar{≈}{\ensuremath{\approx}}\newunicodechar{∂}{\ensuremath{\partial}}\newunicodechar{⋈}{\ensuremath{\bowtie}}\newunicodechar{Δ}{\ensuremath{\Delta}}\newunicodechar{τ}{\ensuremath{\tau}}\newunicodechar{σ}{\ensuremath{\sigma}}\newunicodechar{β}{\ensuremath{\beta}}\newunicodechar{θ}{\ensuremath{\theta}}\newunicodechar{γ}{\ensuremath{\gamma}}\newunicodechar{α}{\ensuremath{\alpha}}\newunicodechar{δ}{\ensuremath{\delta}}\newunicodechar{ε}{\ensuremath{\varepsilon}}\newunicodechar{²}{\textsuperscript{2}}\newunicodechar{⇒}{\ensuremath{\Rightarrow}}\newunicodechar{▶}{\ensuremath{\blacktriangleright}}\newunicodechar{│}{|}\newunicodechar{▼}{\ensuremath{\blacktriangledown}}\newunicodechar{─}{-}\newunicodechar{┌}{+}\newunicodechar{┐}{+}\newunicodechar{└}{+}\newunicodechar{┘}{+}\newunicodechar{├}{+}\newunicodechar{┤}{+} | |
| 8 | +\usepackage{xcolor} | |
| 9 | +\definecolor{uqo}{HTML}{006085}\definecolor{turq}{HTML}{2BCCAE}\definecolor{corail}{HTML}{EA6047}\definecolor{encre}{HTML}{16181B}\definecolor{papier}{HTML}{F4F5F3} | |
| 10 | +\usepackage{titlesec} | |
| 11 | +\titleformat{\section}{\Large\bfseries\color{uqo}}{}{0em}{}[{\color{uqo}\titlerule[0.8pt]}] | |
| 12 | +\titleformat{\subsection}{\large\bfseries\color{uqo}}{}{0em}{} | |
| 13 | +\titlespacing*{\section}{0pt}{18pt plus 4pt}{8pt} | |
| 14 | +\usepackage{etoolbox} | |
| 15 | +\AtBeginEnvironment{longtable}{\scriptsize\renewcommand{\arraystretch}{1.18}} | |
| 16 | +\usepackage{colortbl}\arrayrulecolor{uqo!35} | |
| 17 | +\usepackage{tcolorbox} | |
| 18 | +\renewenvironment{quote}{\begin{tcolorbox}[colback=papier,colframe=turq,boxrule=0pt,leftrule=3pt,arc=1pt,left=6pt,right=6pt,top=4pt,bottom=4pt]}{\end{tcolorbox}} | |
| 19 | +\usepackage{fancyhdr} | |
| 20 | +\pagestyle{fancy}\fancyhf{} | |
| 21 | +\fancyhead[L]{\footnotesize\color{gray}Grille empirique d'ajustement des terrains — UQO Éval}\fancyhead[R]{\footnotesize\color{gray}IMM1003 · IMM1033} | |
| 22 | +\fancyfoot[C]{\footnotesize\color{gray}\thepage}\renewcommand{\headrulewidth}{0.3pt} | |
| 23 | +\usepackage{caption}\captionsetup{font=small,labelfont={bf,color=uqo}} | |
| 24 | +\usepackage{float}\floatplacement{figure}{H} | |
| 25 | +\setlength{\emergencystretch}{3em}\hyphenpenalty=300\tolerance=1500 | |
| 26 | +\usepackage{microtype} | |
| 27 | +\makeatletter\renewcommand{\maketitle}{\begin{center}{\color{uqo}\Huge\bfseries\@title\par}\vspace{8pt}{\normalsize\color{gray}\@author\par\@date\par}\end{center}\vspace{10pt}}\makeatother | |
| 28 | +\usepackage{polyglossia}\setdefaultlanguage[frenchpart=false]{french} | |
added
make_pdf.sh
+13 −0
@@ -0,0 +1,13 @@ | ||
| 1 | +#!/bin/zsh | |
| 2 | +# Rapatrie les résultats depuis M2U64 et produit docs/RAPPORT.pdf (pandoc → XeLaTeX) | |
| 3 | +set -e | |
| 4 | +cd "$(dirname "$0")" | |
| 5 | +V=${CV_VERSION:-2026-10-10} | |
| 6 | +mkdir -p out/$V | |
| 7 | +rsync -az --exclude 'est_E*.npy' --exclude 'sd_E2.npy' --exclude '*_pred_*.npy' --exclude 'sales.parquet' M2U64:comparables-validation/out/$V/ out/$V/ | |
| 8 | +rsync -az M2U64:comparables-validation/docs/ docs/ 2>/dev/null || true | |
| 9 | +cd docs | |
| 10 | +pandoc RAPPORT.md -o RAPPORT.pdf --pdf-engine=xelatex -H header.tex --toc --toc-depth=2 \ | |
| 11 | + -V geometry:"a4paper,margin=17mm,top=20mm,bottom=20mm" -V fontsize=10pt -V linkcolor=uqo -V urlcolor=uqo -V toccolor=uqo -V colorlinks=true \ | |
| 12 | + -V toc-title="Table des matières" --resource-path=.:../out/$V/fig 2>&1 | grep -v "Missing character\|Overfull\|Underfull" | tail -8 || true | |
| 13 | +echo "PDF : $(pwd)/RAPPORT.pdf ($(du -h RAPPORT.pdf | cut -f1))" | |
added
pipeline/01_dataset.py
+97 −0
@@ -0,0 +1,97 @@ | ||
| 1 | +# -*- coding: utf-8 -*- | |
| 2 | +"""01 — Jeu de données : transactions unifamiliales Vrai Prix ⋈ rôle, nettoyage journalisé, | |
| 3 | +reventes, géographie, périodes. Sorties : sales.parquet, exclusions.csv, dataset_stats.json.""" | |
| 4 | +from common import * | |
| 5 | + | |
| 6 | +con = connect() | |
| 7 | +log(f"lecture {DB}") | |
| 8 | +sql = """ | |
| 9 | +select t.id, t.date, t.amount, t.lat, t.lng, t.city, t.property_type, t.year_built, t.floor_area, | |
| 10 | + t.building_type, t.id_provinc, t.valeur_role as role_vente, t.land_area, | |
| 11 | + u.code_mun, u.municipalite, u.type_prop as u_type, u.annee_construction, u.aire_etages_m2, | |
| 12 | + u.superficie_terrain_m2, u.front_terrain_m, u.nb_etages, u.nb_logements, u.lien_physique, | |
| 13 | + u.genre_construction, u.unite_voisinage, u.dat_cond_marche, u.valeur_terrain, u.valeur_batiment, | |
| 14 | + u.valeur_role as role_actuel, u.valeur_anterieure as role_anterieur, u.cubf | |
| 15 | +from transactions t left join units u on u.id_provinc = t.id_provinc | |
| 16 | +where t.property_type = 'unifamilial' | |
| 17 | +""" | |
| 18 | +df = pd.read_sql(sql, con) | |
| 19 | +log(f"{len(df):,} transactions unifamiliales brutes") | |
| 20 | +n0 = len(df) | |
| 21 | +excl = [] # (règle, n) | |
| 22 | +def drop(mask, rule): | |
| 23 | + global df | |
| 24 | + n = int(mask.sum()); excl.append((rule, n)); df = df[~mask].copy(); log(f" {rule}: -{n:,} → {len(df):,}") | |
| 25 | + | |
| 26 | +df["date"] = pd.to_datetime(df["date"]) | |
| 27 | +df = df.sort_values(["id_provinc", "date"]).reset_index(drop=True) | |
| 28 | +drop(df["id_provinc"].isna() | df["code_mun"].isna(), "R0 sans jumelage au rôle (id_provinc / code_mun)") | |
| 29 | +drop(df.duplicated(["id_provinc", "date"]), "R1 doublons (même unité, même date)") | |
| 30 | +drop(~df["u_type"].isin(["unifamilial"]) | (df["nb_logements"].fillna(1) > 1), "R2 unité du rôle non unifamiliale (chalet, plex, terrain…) ou > 1 logement") | |
| 31 | +drop(df["lat"].isna() | df["lng"].isna() | (df["lat"] < 44.5) | (df["lat"] > 63) | (df["lng"] > -56) | (df["lng"] < -80), "R3 coordonnées absentes ou hors Québec") | |
| 32 | +fa = df["floor_area"].fillna(df["aire_etages_m2"]) | |
| 33 | +df["floor_area"] = fa | |
| 34 | +drop(fa.isna() | (fa < 40) | (fa > 600), "R4 aire habitable absente ou hors [40, 600] m²") | |
| 35 | +la = df["land_area"].fillna(df["superficie_terrain_m2"]) | |
| 36 | +df["land_area"] = la | |
| 37 | +drop(la.isna() | (la < 80) | (la > 20000), "R5 superficie de terrain absente ou hors [80, 20 000] m²") | |
| 38 | +yb = df["year_built"].fillna(df["annee_construction"]) | |
| 39 | +df["year_built"] = yb | |
| 40 | +drop(yb.isna() | (yb < 1850) | (yb > df["date"].dt.year), "R6 année de construction absente, < 1850 ou postérieure à la vente (vente de terrain)") | |
| 41 | +drop((df["year_built"] == df["date"].dt.year), "R6b vente l'année même de la construction (prix du terrain ou d'une construction inachevée : médiane 346 k$ et ratio au rôle 1,15 contre 450-480 k$ et 1,3 aux âges 1-4)") | |
| 42 | +drop((df["amount"] < 50000) | (df["amount"] > 3_000_000), "R7 prix hors [50 k$, 3 M$]") | |
| 43 | +ratio = df["amount"] / df["role_vente"] | |
| 44 | +drop(df["role_vente"].isna() | (df["role_vente"] <= 0) | (ratio < 0.33) | (ratio > 3.0), "R8 rôle absent ou ratio prix/rôle hors [0,33 ; 3] (transactions présumées non représentatives)") | |
| 45 | +# reventes rapides avec forte variation (flips / rénovations : les caractéristiques enregistrées ne décrivent plus le bien) | |
| 46 | +df = df.sort_values(["id_provinc", "date"]).reset_index(drop=True) | |
| 47 | +g = df.groupby("id_provinc") | |
| 48 | +prev_date = g["date"].shift(1); prev_amt = g["amount"].shift(1) | |
| 49 | +dd = (df["date"] - prev_date).dt.days | |
| 50 | +chg = df["amount"] / prev_amt - 1 | |
| 51 | +flip = (dd <= 180) & (chg.abs() > 0.40) | |
| 52 | +drop(flip.fillna(False), "R9 revente ≤ 180 jours avec variation > 40 % (flip)") | |
| 53 | +drop(df["lien_physique"].isna() | df["genre_construction"].isna(), "R10 lien physique ou genre de construction absent au rôle") | |
| 54 | + | |
| 55 | +# ----------------------------------------------------------------- variables dérivées | |
| 56 | +rep = load_repertoire() | |
| 57 | +df = df.merge(rep, on="code_mun", how="left") | |
| 58 | +drop(df["region"].isna(), "R11 municipalité absente du répertoire MAMH") | |
| 59 | +df["sale_year"] = df["date"].dt.year | |
| 60 | +df["month"] = month_index(df["date"]) | |
| 61 | +df["day"] = (df["date"] - pd.Timestamp("2021-01-01")).dt.days.astype(int) | |
| 62 | +df["age"] = (df["sale_year"] - df["year_built"]).clip(lower=0) | |
| 63 | +df["ln_price"] = np.log(df["amount"]) | |
| 64 | +df["ln_floor"] = np.log(df["floor_area"]) | |
| 65 | +df["ln_land"] = np.log(df["land_area"]) | |
| 66 | +df["lien"] = df["lien_physique"].map(LIEN).fillna(0).astype(int) | |
| 67 | +df["genre"] = df["genre_construction"].map(GENRE).fillna(0).astype(int) | |
| 68 | +df["btype"] = df["building_type"].map(BTYPE).fillna(-1).astype(int) | |
| 69 | +df["etages"] = df["nb_etages"].fillna(1).clip(1, 3) | |
| 70 | +df["x_km"], df["y_km"] = to_xy_km(df["lat"], df["lng"]) | |
| 71 | +# millésime du rôle au moment de la vente (rôle courant / antérieur / autre) | |
| 72 | +df["role_vintage"] = np.where(np.isclose(df["role_vente"], df["role_actuel"].fillna(-1)), 2, | |
| 73 | + np.where(np.isclose(df["role_vente"], df["role_anterieur"].fillna(-1)), 1, 0)) | |
| 74 | +df["ratio_role"] = df["amount"] / df["role_vente"] | |
| 75 | +df["part_terrain"] = (df["valeur_terrain"] / df["role_actuel"]).clip(0.02, 0.95).fillna(0.2) | |
| 76 | +df["neuf"] = (df["age"] <= 3).astype(int) | |
| 77 | +df["period"] = np.select([df["date"] <= TRAIN_END, df["date"] <= VAL_END], ["train", "val"], "test") | |
| 78 | +df = df.sort_values("date").reset_index(drop=True) | |
| 79 | +df["idx"] = np.arange(len(df), dtype=np.int32) | |
| 80 | + | |
| 81 | +keep = ["idx", "id", "id_provinc", "date", "day", "month", "sale_year", "period", "amount", "ln_price", | |
| 82 | + "lat", "lng", "x_km", "y_km", "code_mun", "mun_nom", "mrc", "mrc_nom", "region", "region_nom", | |
| 83 | + "unite_voisinage", "floor_area", "ln_floor", "land_area", "ln_land", "year_built", "age", "neuf", | |
| 84 | + "lien", "genre", "btype", "etages", "role_vente", "role_vintage", "ratio_role", "part_terrain"] | |
| 85 | +df = df[keep] | |
| 86 | +save_parquet(df, "sales.parquet") | |
| 87 | +pd.DataFrame(excl, columns=["regle", "n"]).to_csv(OUT / "exclusions.csv", index=False) | |
| 88 | +stats = { | |
| 89 | + "brut": n0, "propre": len(df), "exclusions": excl, | |
| 90 | + "par_periode": df["period"].value_counts().to_dict(), | |
| 91 | + "par_annee": df["sale_year"].value_counts().sort_index().to_dict(), | |
| 92 | + "par_region": df.groupby("region_nom").size().sort_values(ascending=False).to_dict(), | |
| 93 | + "prix": {k: float(v) for k, v in df["amount"].describe().items()}, | |
| 94 | + "dates": [str(df["date"].min().date()), str(df["date"].max().date())], | |
| 95 | +} | |
| 96 | +save_json(stats, "dataset_stats.json") | |
| 97 | +log(f"jeu propre : {len(df):,} ventes ({stats['par_periode']})") | |
added
pipeline/02_index.py
+78 −0
@@ -0,0 +1,78 @@ | ||
| 1 | +# -*- coding: utf-8 -*- | |
| 2 | +"""02 — Indice temporel régional « point dans le temps ». | |
| 3 | +Pour chaque mois d'évaluation m, on ré-estime un indice hédonique mensuel par région administrative | |
| 4 | +avec toutes les ventes antérieures au mois m (fenêtre expansive) (jamais le mois m lui-même). Niveau de référence = m-1, | |
| 5 | +dernier mois entièrement connu. Sortie : index_asof.npy [m_eval, région, mois] (log), regions.json, | |
| 6 | +et un indice « plein échantillon » descriptif (index_full.csv).""" | |
| 7 | +from common import * | |
| 8 | + | |
| 9 | +df = load_parquet("sales.parquet") | |
| 10 | +regions = sorted(df["region"].unique().tolist()) | |
| 11 | +rid = {r: i for i, r in enumerate(regions)} | |
| 12 | +df["rid"] = df["region"].map(rid).astype(int) | |
| 13 | +M = int(df["month"].max()) + 1 | |
| 14 | +M_EVAL0 = int(month_index(pd.Series([VAL_START]))[0]) | |
| 15 | +WINDOW = 999 # fenêtre « expansive » : toutes les ventes antérieures au mois m | |
| 16 | +SHRINK = 40.0 # n/(n+SHRINK) : rétrécissement de l'effet mensuel régional vers l'effet provincial | |
| 17 | + | |
| 18 | +X_cols = ["ln_floor", "ln_land", "age"] | |
| 19 | +def design(d): | |
| 20 | + X = np.column_stack([np.ones(len(d)), d["ln_floor"], d["ln_land"], d["age"], d["age"] ** 2 / 100, | |
| 21 | + (d["lien"] == 1), (d["lien"] >= 2), (d["genre"] == 1), (d["genre"] == 2), (d["genre"] == 3), | |
| 22 | + d["etages"], d["neuf"]]).astype(float) | |
| 23 | + return X | |
| 24 | + | |
| 25 | +def fit_month_effects(d, months): | |
| 26 | + """OLS à effets fixes de municipalité (estimateur within) : ln P ~ X + effets mensuels (base = dernier mois), | |
| 27 | + toutes les variables étant centrées par municipalité pour que l'indice ne soit pas contaminé par les | |
| 28 | + déplacements du mélange des ventes entre municipalités chères et bon marché. Renvoie effets et n par mois.""" | |
| 29 | + X = design(d)[:, 1:] | |
| 30 | + mm = d["month"].to_numpy() | |
| 31 | + base = months[-1] | |
| 32 | + dums = np.column_stack([(mm == m).astype(float) for m in months if m != base]) | |
| 33 | + A = np.hstack([X, dums]); y = d["ln_price"].to_numpy(float) | |
| 34 | + codes, inv = np.unique(d["code_mun"].to_numpy(), return_inverse=True) | |
| 35 | + cnt = np.bincount(inv).astype(float) | |
| 36 | + mA = np.zeros((len(codes), A.shape[1])); np.add.at(mA, inv, A); A = A - mA[inv] / cnt[inv, None] | |
| 37 | + my = np.bincount(inv, y) / cnt; y = y - my[inv] | |
| 38 | + beta, *_ = np.linalg.lstsq(A, y, rcond=None) | |
| 39 | + eff = {base: 0.0} | |
| 40 | + k = X.shape[1] | |
| 41 | + for j, m in enumerate([m for m in months if m != base]): | |
| 42 | + eff[m] = float(beta[k + j]) | |
| 43 | + n = {m: int((mm == m).sum()) for m in months} | |
| 44 | + return eff, n | |
| 45 | + | |
| 46 | +idx = np.full((M - M_EVAL0, len(regions), M), np.nan, dtype=np.float32) | |
| 47 | +log(f"indice point-dans-le-temps : mois d'évaluation {M_EVAL0}..{M-1}, {len(regions)} régions") | |
| 48 | +for m in range(M_EVAL0, M): | |
| 49 | + months = list(range(max(0, m - WINDOW), m)) | |
| 50 | + win = df[(df["month"] >= months[0]) & (df["month"] < m)] | |
| 51 | + prov_eff, _ = fit_month_effects(win, months) | |
| 52 | + for r in regions: | |
| 53 | + dr = win[win["region"] == r] | |
| 54 | + if len(dr) < 200: | |
| 55 | + for mo in months: idx[m - M_EVAL0, rid[r], mo] = prov_eff[mo] | |
| 56 | + continue | |
| 57 | + eff, n = fit_month_effects(dr, months) | |
| 58 | + for mo in months: | |
| 59 | + w = n[mo] / (n[mo] + SHRINK) | |
| 60 | + idx[m - M_EVAL0, rid[r], mo] = w * eff[mo] + (1 - w) * prov_eff[mo] | |
| 61 | + if m % 6 == 0: log(f" mois {m} ({len(win):,} ventes dans la fenêtre)") | |
| 62 | +np.save(OUT / "index_asof.npy", idx) | |
| 63 | +save_json({"regions": regions, "M_EVAL0": M_EVAL0, "M": M, "window": WINDOW, "shrink": SHRINK}, "index_meta.json") | |
| 64 | + | |
| 65 | +# indice descriptif plein échantillon (pour le rapport seulement — pas utilisé dans la rétro-évaluation) | |
| 66 | +rows = [] | |
| 67 | +months = list(range(0, M)) | |
| 68 | +eff_p, n_p = fit_month_effects(df, months) | |
| 69 | +for mo in months: | |
| 70 | + rows.append({"month": mo, "ym": f"{2021 + mo // 12}-{mo % 12 + 1:02d}", "region": "QC", "log_idx": eff_p[mo] - eff_p[months[0]], "n": n_p[mo]}) | |
| 71 | +for r in regions: | |
| 72 | + dr = df[df["region"] == r] | |
| 73 | + if len(dr) < 2000: continue | |
| 74 | + eff, n = fit_month_effects(dr, months) | |
| 75 | + for mo in months: | |
| 76 | + rows.append({"month": mo, "ym": f"{2021 + mo // 12}-{mo % 12 + 1:02d}", "region": r, "log_idx": eff[mo] - eff[months[0]], "n": n[mo]}) | |
| 77 | +pd.DataFrame(rows).to_csv(OUT / "index_full.csv", index=False) | |
| 78 | +log("indice terminé") | |
added
pipeline/03_models.py
+107 −0
@@ -0,0 +1,107 @@ | ||
| 1 | +# -*- coding: utf-8 -*- | |
| 2 | +"""03 — Modèles d'ajustement et modèles de référence, ré-estimés chaque trimestre avec les seules ventes | |
| 3 | +antérieures au trimestre (as-of 2024-01, 2024-04, …, 2026-07) — comme le ferait un AVM ré-entraîné régulièrement : | |
| 4 | + • grille hédonique globale (coefficients) et régionale rétrécie → ajustements B / C | |
| 5 | + • LightGBM sur le prix ajusté dans le temps (avec et sans rôle) → ajustements C-GBM, M11 | |
| 6 | + • hédonique globale directe avec effets de municipalité → M10 | |
| 7 | +Sorties : models_<période>.json, gbm_pred_<période>.npy (prédiction f pour TOUTES les ventes).""" | |
| 8 | +import pickle | |
| 9 | +from common import * | |
| 10 | +import lightgbm as lgb | |
| 11 | + | |
| 12 | +df = load_parquet("sales.parquet") | |
| 13 | +meta = load_json("index_meta.json"); regions = meta["regions"]; M_EVAL0 = meta["M_EVAL0"] | |
| 14 | +rid = {r: i for i, r in enumerate(regions)} | |
| 15 | +df["rid"] = df["region"].map(rid).astype(int) | |
| 16 | +idx = np.load(OUT / "index_asof.npy") | |
| 17 | + | |
| 18 | +FEATS = ["ln_floor", "ln_land", "age", "age2", "lien1", "lien2", "genre1", "genre2", "genre3", "etages", "neuf"] | |
| 19 | +def feats(d): | |
| 20 | + return pd.DataFrame({ | |
| 21 | + "ln_floor": d["ln_floor"], "ln_land": d["ln_land"], "age": d["age"], "age2": d["age"] ** 2 / 100, | |
| 22 | + "lien1": (d["lien"] == 1).astype(float), "lien2": (d["lien"] >= 2).astype(float), | |
| 23 | + "genre1": (d["genre"] == 1).astype(float), "genre2": (d["genre"] == 2).astype(float), "genre3": (d["genre"] == 3).astype(float), | |
| 24 | + "etages": d["etages"].astype(float), "neuf": d["neuf"].astype(float)}, index=d.index) | |
| 25 | + | |
| 26 | +def time_adjusted(d, m_asof): | |
| 27 | + """ln P ramené au niveau du mois m_asof-1 avec l'indice connu au mois m_asof.""" | |
| 28 | + lvl = idx[m_asof - M_EVAL0] | |
| 29 | + return d["ln_price"].to_numpy() - lvl[d["rid"].to_numpy(), d["month"].to_numpy()] + lvl[d["rid"].to_numpy(), m_asof - 1] | |
| 30 | + | |
| 31 | +def ols(X, y): | |
| 32 | + A = np.column_stack([np.ones(len(X)), X]) | |
| 33 | + b, *_ = np.linalg.lstsq(A, y, rcond=None) | |
| 34 | + return b | |
| 35 | + | |
| 36 | +ASOF = list(range(M_EVAL0, int(df["month"].max()) + 1, 3)) # ré-estimation trimestrielle (2024-01, 2024-04, …) | |
| 37 | +save_json({"asof": ASOF}, "asof.json") | |
| 38 | +for m_asof in ASOF: | |
| 39 | + period = str(m_asof); start = f"{2021 + m_asof // 12}-{m_asof % 12 + 1:02d}-01" | |
| 40 | + tr = df[df["month"] < m_asof].copy() | |
| 41 | + y = time_adjusted(tr, m_asof) | |
| 42 | + F = feats(tr) | |
| 43 | + log(f"[{period}] entraînement sur {len(tr):,} ventes < {start} (indice connu au mois {m_asof})") | |
| 44 | + # --- grille hédonique globale : estimateur « within » (effets fixes de municipalité), pour que les β | |
| 45 | + # mesurent des différences entre biens voisins et non entre marchés | |
| 46 | + def within(frame, cols, by): | |
| 47 | + return frame[cols] - frame.groupby(by)[cols].transform("mean") | |
| 48 | + Fw = within(pd.concat([F, tr[["code_mun"]]], axis=1), FEATS, "code_mun") | |
| 49 | + yw = y - pd.Series(y, index=tr.index).groupby(tr["code_mun"]).transform("mean").to_numpy() | |
| 50 | + bw, *_ = np.linalg.lstsq(Fw.to_numpy(), yw, rcond=None) | |
| 51 | + beta_global = dict(zip(FEATS, bw.tolist())) | |
| 52 | + xb = F.to_numpy() @ bw | |
| 53 | + resid = y - xb | |
| 54 | + r2w = 1 - np.var(yw - Fw.to_numpy() @ bw) / np.var(yw) | |
| 55 | + log(f" hédonique within-municipalité : R² within = {r2w:.3f} ; β = { {k: round(v, 3) for k, v in beta_global.items()} }") | |
| 56 | + # --- grille régionale (within-municipalité) rétrécie vers la globale (n/(n+300)) | |
| 57 | + beta_region = {} | |
| 58 | + for r in regions: | |
| 59 | + sel = (tr["region"] == r).to_numpy() | |
| 60 | + if sel.sum() < 100: beta_region[r] = beta_global; continue | |
| 61 | + b, *_ = np.linalg.lstsq(Fw.to_numpy()[sel], yw[sel], rcond=None) | |
| 62 | + w = sel.sum() / (sel.sum() + 300.0) | |
| 63 | + beta_region[r] = {k: float(w * b[j] + (1 - w) * beta_global[k]) for j, k in enumerate(FEATS)} | |
| 64 | + # --- M10 : hédonique directe = Xβ + effet de région + MRC + municipalité + unité de voisinage (rétrécis en cascade) | |
| 65 | + reg_eff = pd.Series(resid, index=tr.index).groupby(tr["region"]).mean().to_dict() | |
| 66 | + tr["res_g"] = resid - tr["region"].map(reg_eff).to_numpy() | |
| 67 | + mrc_eff = tr.groupby("mrc")["res_g"].agg(["sum", "count"]) | |
| 68 | + mrc_eff = (mrc_eff["sum"] / (mrc_eff["count"] + 30)).to_dict() | |
| 69 | + tr["res_m"] = tr["res_g"] - tr["mrc"].map(mrc_eff).fillna(0) | |
| 70 | + mun_eff = tr.groupby("code_mun")["res_m"].agg(["sum", "count"]) | |
| 71 | + mun_eff = (mun_eff["sum"] / (mun_eff["count"] + 20)).to_dict() | |
| 72 | + tr["res_u"] = tr["res_m"] - tr["code_mun"].map(mun_eff).fillna(0) | |
| 73 | + uv = tr.groupby(["code_mun", "unite_voisinage"])["res_u"].agg(["sum", "count"]) | |
| 74 | + uv_eff = {f"{a}|{b}": v for (a, b), v in (uv["sum"] / (uv["count"] + 10)).items()} | |
| 75 | + # --- LightGBM (prix ajusté dans le temps) avec rôle et sans rôle | |
| 76 | + def gbm_frame(d): | |
| 77 | + G = feats(d) | |
| 78 | + G["lat"] = d["lat"]; G["lng"] = d["lng"]; G["x_km"] = d["x_km"]; G["y_km"] = d["y_km"] | |
| 79 | + G["rid"] = d["rid"]; G["btype"] = d["btype"] | |
| 80 | + G["ln_role"] = np.log(d["role_vente"]); G["part_terrain"] = d["part_terrain"]; G["role_vintage"] = d["role_vintage"] | |
| 81 | + return G | |
| 82 | + G = gbm_frame(tr) | |
| 83 | + params = dict(objective="regression", learning_rate=0.05, num_leaves=127, min_data_in_leaf=40, feature_fraction=0.8, | |
| 84 | + bagging_fraction=0.8, bagging_freq=1, lambda_l2=1.0, verbose=-1, num_threads=20) | |
| 85 | + models = {} | |
| 86 | + for name, cols in [("gbm", [c for c in G.columns if c not in ("ln_role", "part_terrain", "role_vintage")]), | |
| 87 | + ("gbm_role", list(G.columns))]: | |
| 88 | + ds = lgb.Dataset(G[cols], y, categorical_feature=["rid", "btype"] if "btype" in cols else ["rid"]) | |
| 89 | + # validation interne : derniers 10 % chronologiques pour fixer le nombre d'arbres | |
| 90 | + n = len(G); cut = int(n * 0.9) | |
| 91 | + dtr = lgb.Dataset(G[cols].iloc[:cut], y[:cut], categorical_feature=["rid", "btype"]) | |
| 92 | + dva = lgb.Dataset(G[cols].iloc[cut:], y[cut:], reference=dtr) | |
| 93 | + bst = lgb.train(params, dtr, num_boost_round=3000, valid_sets=[dva], callbacks=[lgb.early_stopping(100, verbose=False)]) | |
| 94 | + best = bst.best_iteration | |
| 95 | + bst = lgb.train(params, ds, num_boost_round=best) | |
| 96 | + models[name] = bst | |
| 97 | + pred_all = bst.predict(gbm_frame(df)[cols]) | |
| 98 | + np.save(OUT / f"{name}_pred_{period}.npy", pred_all.astype(np.float32)) | |
| 99 | + r2 = 1 - np.var(y - bst.predict(G[cols])) / np.var(y) | |
| 100 | + log(f" {name}: {best} arbres, R² (in-sample) = {r2:.3f}") | |
| 101 | + bst.save_model(str(OUT / f"{name}_{period}.txt")) | |
| 102 | + save_json({"period": period, "m_asof": m_asof, "start": start, "n_train": len(tr), "feats": FEATS, | |
| 103 | + "beta_global": beta_global, "intercept_global": 0.0, "region_fe": reg_eff, | |
| 104 | + "beta_region": beta_region, "mrc_eff": mrc_eff, "mun_eff": mun_eff, "uv_eff": uv_eff, | |
| 105 | + "r2_within": float(r2w), "sd_resid_u": float(tr["res_u"].std())}, | |
| 106 | + f"models_{period}.json") | |
| 107 | +log("modèles terminés") | |
added
pipeline/04_backtest.py
+190 −0
@@ -0,0 +1,190 @@ | ||
| 1 | +# -*- coding: utf-8 -*- | |
| 2 | +"""04 — Moteur de rétro-évaluation. | |
| 3 | +Pour chaque vente cible (2024 = validation, 2025-2026 = test), le prix est masqué ; seules les ventes | |
| 4 | +conclues STRICTEMENT AVANT la date cible (et l'indice/les modèles connus à cette date) sont utilisées. | |
| 5 | +On évalue d'un coup toutes les configurations (sélection × rayon × fenêtre × nombre × ajustement × | |
| 6 | +pondération × réconciliation) sur exactement les mêmes cibles. | |
| 7 | +Sorties : targets.parquet, configs_E1.csv/est_E1.npy (surface k,d,T), configs_E2.csv/est_E2.npy/sd_E2.npy | |
| 8 | +(factoriel sélection × ajustement × pondération × réconciliation), n_avail.npy.""" | |
| 9 | +import os | |
| 10 | +for _v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"): | |
| 11 | + os.environ[_v] = "1" # un seul fil BLAS par processus (fork + BLAS multi-fils = interblocage sur macOS) | |
| 12 | +import multiprocessing as mp, time | |
| 13 | +from scipy.spatial import cKDTree | |
| 14 | +from common import * | |
| 15 | + | |
| 16 | +df = load_parquet("sales.parquet") | |
| 17 | +meta = load_json("index_meta.json"); regions = meta["regions"]; M_EVAL0 = meta["M_EVAL0"] | |
| 18 | +rid_map = {r: i for i, r in enumerate(regions)} | |
| 19 | +IDX = np.load(OUT / "index_asof.npy") | |
| 20 | +ASOF = load_json("asof.json")["asof"] | |
| 21 | +MODELS = {a: load_json(f"models_{a}.json") for a in ASOF} | |
| 22 | +FEATS = MODELS[ASOF[0]]["feats"] | |
| 23 | +GBM = {a: np.load(OUT / f"gbm_pred_{a}.npy") for a in ASOF} | |
| 24 | + | |
| 25 | +# ----------------------------------------------------------------- tableaux globaux (partagés par fork) | |
| 26 | +N = len(df) | |
| 27 | +DAY = df["day"].to_numpy(np.int32); MONTH = df["month"].to_numpy(np.int32) | |
| 28 | +RID = df["region"].map(rid_map).to_numpy(np.int32) | |
| 29 | +X = df["x_km"].to_numpy(np.float64); Y = df["y_km"].to_numpy(np.float64) | |
| 30 | +LNP = df["ln_price"].to_numpy(np.float64); P = df["amount"].to_numpy(np.float64) | |
| 31 | +FLOOR = df["floor_area"].to_numpy(np.float64); LAND = df["land_area"].to_numpy(np.float64) | |
| 32 | +AGE = df["age"].to_numpy(np.float64); LIEN = df["lien"].to_numpy(np.int8); ETAGES = df["etages"].to_numpy(np.float64) | |
| 33 | +GENRE = df["genre"].to_numpy(np.int8) | |
| 34 | +PROV = pd.factorize(df["id_provinc"])[0].astype(np.int32) | |
| 35 | +MUN = pd.factorize(df["code_mun"])[0].astype(np.int32) | |
| 36 | +VINT = df["role_vintage"].to_numpy(np.int8); LNROLE = np.log(df["role_vente"].to_numpy(np.float64)) | |
| 37 | +PART = df["part_terrain"].to_numpy(np.float64) | |
| 38 | +NEUF = df["neuf"].to_numpy(np.float64) | |
| 39 | +# matrice des caractéristiques hédoniques (même ordre que FEATS) | |
| 40 | +XF = np.column_stack([np.log(FLOOR), np.log(LAND), AGE, AGE ** 2 / 100, (LIEN == 1), (LIEN >= 2), | |
| 41 | + (GENRE == 1), (GENRE == 2), (GENRE == 3), ETAGES, NEUF]).astype(np.float64) | |
| 42 | +BETA_G = {a: np.array([MODELS[a]["beta_global"][f] for f in FEATS]) for a in ASOF} | |
| 43 | +BETA_R = {a: np.array([[MODELS[a]["beta_region"][r][f] for f in FEATS] for r in regions]) for a in ASOF} | |
| 44 | +ASOF_ARR = np.array(ASOF) | |
| 45 | +TREE = cKDTree(np.column_stack([X, Y])) | |
| 46 | + | |
| 47 | +D_MAX, T_MAX_DAYS = 25.0, 731 | |
| 48 | +KS = [1, 3, 5, 7, 10, 15, 20, 30, 50] | |
| 49 | +DS = [0.5, 1, 2, 5, 10, 25] | |
| 50 | +TS = [3, 6, 12, 24] | |
| 51 | +SELS = ["geo", "sim", "mix", "adj"] | |
| 52 | +ADJS = ["raw", "none", "conv", "hedo", "hedo_reg", "local", "gbm", "role"] | |
| 53 | +WRS = [("equal", "median"), ("equal", "mean"), ("dist", "median"), ("dist", "mean"), ("rec", "median"), ("rec", "mean"), | |
| 54 | + ("sim", "median"), ("sim", "mean"), ("all", "median"), ("all", "mean")] | |
| 55 | +ANCHORS = [(3, 2, 6), (5, 5, 12), (10, 10, 24)] | |
| 56 | +E1 = [dict(sel=s, d=d, T=T, k=k, adj=a, w="equal", recon="median") | |
| 57 | + for s in ("geo", "sim", "mix") for d in DS for T in TS for k in KS for a in ("none", "hedo")] | |
| 58 | +E2 = [dict(sel=s, d=d, T=T, k=k, adj=a, w=w, recon=r) | |
| 59 | + for (k, d, T) in ANCHORS for s in SELS for a in ADJS for (w, r) in WRS] | |
| 60 | +LOCAL_COLS = [0, 1, 2, 4, 5, 9] # ln_floor, ln_land, age, lien1, lien2, etages | |
| 61 | +LAMBDA_LOCAL = 50.0 | |
| 62 | + | |
| 63 | +def wmedian(v, w): | |
| 64 | + o = np.argsort(v); v = v[o]; w = w[o] | |
| 65 | + c = np.cumsum(w); return v[np.searchsorted(c, 0.5 * c[-1])] | |
| 66 | + | |
| 67 | +def recon(vals, w, how): | |
| 68 | + ok = ~np.isnan(vals) | |
| 69 | + if not ok.any(): return np.nan, np.nan | |
| 70 | + v = vals[ok]; ww = w[ok] | |
| 71 | + sd = v.std() if len(v) > 1 else np.nan | |
| 72 | + if how == "median": | |
| 73 | + return (np.median(v) if ww.min() == ww.max() else wmedian(v, ww)), sd | |
| 74 | + return np.log(np.sum(ww * np.exp(v)) / ww.sum()), sd | |
| 75 | + | |
| 76 | +def evaluate(i): | |
| 77 | + """Toutes les configurations pour la cible i. Renvoie (est_E1, est_E2, sd_E2, n_avail, info).""" | |
| 78 | + period = int(ASOF_ARR[ASOF_ARR <= MONTH[i]].max()) # dernier trimestre de ré-estimation ≤ mois de la cible | |
| 79 | + e1 = np.full(len(E1), np.nan, np.float32); e2 = np.full(len(E2), np.nan, np.float32) | |
| 80 | + s2 = np.full(len(E2), np.nan, np.float32); nav = np.zeros(len(DS) * len(TS), np.int32) | |
| 81 | + cand = np.asarray(TREE.query_ball_point((X[i], Y[i]), D_MAX), dtype=np.int64) | |
| 82 | + keep = (DAY[cand] < DAY[i]) & (DAY[cand] >= DAY[i] - T_MAX_DAYS) & (PROV[cand] != PROV[i]) | |
| 83 | + cand = cand[keep] | |
| 84 | + info = dict(idx=int(i), asof=period, n_cand=int(len(cand))) | |
| 85 | + if len(cand) == 0: return e1, e2, s2, nav, info | |
| 86 | + dist = np.hypot(X[cand] - X[i], Y[cand] - Y[i]) | |
| 87 | + mago = (DAY[i] - DAY[cand]) / 30.4375 | |
| 88 | + m_eval = MONTH[i]; lvl = IDX[m_eval - M_EVAL0] | |
| 89 | + tadj = np.where(MONTH[cand] == m_eval, 0.0, lvl[RID[cand], m_eval - 1] - lvl[RID[cand], MONTH[cand]]) | |
| 90 | + tadj = np.nan_to_num(tadj, nan=0.0) | |
| 91 | + lnT = LNP[cand] + tadj | |
| 92 | + dX = XF[i] - XF[cand] # cible − comparable | |
| 93 | + # --- ajustements | |
| 94 | + bg = BETA_G[period]; br = BETA_R[period][RID[i]] | |
| 95 | + adj = {} | |
| 96 | + adj["raw"] = LNP[cand] | |
| 97 | + adj["none"] = lnT | |
| 98 | + Pt = np.exp(lnT) | |
| 99 | + d_floor = (FLOOR[i] - FLOOR[cand]) * 0.5 * Pt / FLOOR[cand] | |
| 100 | + d_land = np.clip((LAND[i] - LAND[cand]) * 0.5 * Pt * PART[cand] / LAND[cand], -0.15 * Pt, 0.15 * Pt) | |
| 101 | + d_age = np.clip(-0.005 * Pt * (AGE[i] - AGE[cand]), -0.2 * Pt, 0.2 * Pt) | |
| 102 | + d_lien = np.where((LIEN[i] == 0) & (LIEN[cand] > 0), 0.05 * Pt, np.where((LIEN[i] > 0) & (LIEN[cand] == 0), -0.05 * Pt, 0.0)) | |
| 103 | + tot = np.maximum(d_floor + d_land + d_age + d_lien, -0.6 * Pt) | |
| 104 | + adj["conv"] = np.log(Pt + tot) | |
| 105 | + gross_conv = (np.abs(d_floor) + np.abs(d_land) + np.abs(d_age) + np.abs(d_lien)) / Pt | |
| 106 | + adj["hedo"] = lnT + dX @ bg | |
| 107 | + gross_hedo = np.abs(dX * bg).sum(axis=1) | |
| 108 | + adj["hedo_reg"] = lnT + dX @ br | |
| 109 | + loc = dist <= 10.0 | |
| 110 | + if loc.sum() >= 40: | |
| 111 | + Xl = XF[cand][loc][:, LOCAL_COLS]; yl = lnT[loc] | |
| 112 | + b0 = br[LOCAL_COLS] | |
| 113 | + Xc = Xl - Xl.mean(axis=0); yc = yl - yl.mean() | |
| 114 | + A = Xc.T @ Xc + LAMBDA_LOCAL * np.eye(len(LOCAL_COLS)) | |
| 115 | + bl = np.linalg.solve(A, Xc.T @ yc + LAMBDA_LOCAL * b0) | |
| 116 | + adj["local"] = lnT + dX[:, LOCAL_COLS] @ bl | |
| 117 | + else: | |
| 118 | + adj["local"] = adj["hedo_reg"] | |
| 119 | + f = GBM[period] | |
| 120 | + adj["gbm"] = lnT + (f[i] - f[cand]) | |
| 121 | + same = (MUN[cand] == MUN[i]) & (VINT[cand] == VINT[i]) | |
| 122 | + adj["role"] = np.where(same, LNROLE[i] + (lnT - LNROLE[cand]), np.nan) | |
| 123 | + # --- scores de sélection | |
| 124 | + dissim = np.sqrt((dX[:, 0] / 0.25) ** 2 + (dX[:, 1] / 0.6) ** 2 + (dX[:, 2] / 15) ** 2 | |
| 125 | + + (LIEN[cand] != LIEN[i]) * 1.0 + (ETAGES[cand] != ETAGES[i]) * 0.5) | |
| 126 | + score = {"geo": dist, "sim": dissim + dist / 25.0, "mix": dissim + dist / 2.0, "adj": gross_hedo + dist / 25.0} | |
| 127 | + # --- pondérations | |
| 128 | + wts = {"equal": np.ones(len(cand)), "dist": 1.0 / (0.25 + dist), "rec": np.exp(-mago / 6.0), "sim": 1.0 / (0.1 + gross_hedo)} | |
| 129 | + wts["all"] = wts["dist"] * wts["rec"] * wts["sim"] | |
| 130 | + # --- masques (d, T) et classements | |
| 131 | + top = {} | |
| 132 | + for di, d in enumerate(DS): | |
| 133 | + for ti, T in enumerate(TS): | |
| 134 | + m = (dist <= d) & (mago <= T) | |
| 135 | + nav[di * len(TS) + ti] = int(m.sum()) | |
| 136 | + if not m.any(): continue | |
| 137 | + ids = np.flatnonzero(m) | |
| 138 | + for s in SELS: | |
| 139 | + sc = score[s][ids] | |
| 140 | + if len(ids) > 50: | |
| 141 | + part = np.argpartition(sc, 50)[:50]; part = part[np.argsort(sc[part])] | |
| 142 | + else: | |
| 143 | + part = np.argsort(sc) | |
| 144 | + top[(s, d, T)] = ids[part] | |
| 145 | + # --- E1 | |
| 146 | + for c, cfg in enumerate(E1): | |
| 147 | + t = top.get((cfg["sel"], cfg["d"], cfg["T"])) | |
| 148 | + if t is None: continue | |
| 149 | + sub = t[:cfg["k"]] | |
| 150 | + e1[c], _ = recon(adj[cfg["adj"]][sub], wts["equal"][sub], "median") | |
| 151 | + # --- E2 | |
| 152 | + for c, cfg in enumerate(E2): | |
| 153 | + t = top.get((cfg["sel"], cfg["d"], cfg["T"])) | |
| 154 | + if t is None: continue | |
| 155 | + sub = t[:cfg["k"]] | |
| 156 | + e2[c], s2[c] = recon(adj[cfg["adj"]][sub], wts[cfg["w"]][sub], cfg["recon"]) | |
| 157 | + return e1, e2, s2, nav, info | |
| 158 | + | |
| 159 | +def worker(chunk): | |
| 160 | + return [evaluate(i) for i in chunk] | |
| 161 | + | |
| 162 | +if __name__ == "__main__": | |
| 163 | + targets = df.index[df["period"].isin(["val", "test"])].to_numpy() | |
| 164 | + if os.environ.get("CV_SAMPLE"): | |
| 165 | + rng = np.random.default_rng(1); targets = np.sort(rng.choice(targets, int(os.environ["CV_SAMPLE"]), replace=False)) | |
| 166 | + log(f"{len(targets):,} cibles ; {len(E1)} configurations E1, {len(E2)} E2 ; {N:,} ventes dans le bassin") | |
| 167 | + pd.DataFrame(E1).to_csv(OUT / "configs_E1.csv", index_label="cfg") | |
| 168 | + pd.DataFrame(E2).to_csv(OUT / "configs_E2.csv", index_label="cfg") | |
| 169 | + est1 = np.full((len(targets), len(E1)), np.nan, np.float32) | |
| 170 | + est2 = np.full((len(targets), len(E2)), np.nan, np.float32) | |
| 171 | + sd2 = np.full((len(targets), len(E2)), np.nan, np.float32) | |
| 172 | + nav = np.zeros((len(targets), len(DS) * len(TS)), np.int32) | |
| 173 | + infos = [] | |
| 174 | + chunks = [targets[j:j + 400] for j in range(0, len(targets), 400)] | |
| 175 | + t0 = time.time(); done = 0 | |
| 176 | + with mp.get_context(os.environ.get("CV_MP", "fork")).Pool(int(os.environ.get("CV_PROCS", 20))) as pool: | |
| 177 | + for res in pool.imap(worker, chunks): | |
| 178 | + for (e1, e2, s2, nv, info) in res: | |
| 179 | + est1[done] = e1; est2[done] = e2; sd2[done] = s2; nav[done] = nv; infos.append(info); done += 1 | |
| 180 | + if done % 20000 < 400: | |
| 181 | + log(f" {done:,}/{len(targets):,} cibles ({time.time() - t0:.0f} s)") | |
| 182 | + np.save(OUT / "est_E1.npy", est1); np.save(OUT / "est_E2.npy", est2); np.save(OUT / "sd_E2.npy", sd2); np.save(OUT / "n_avail.npy", nav) | |
| 183 | + T = df.loc[targets, ["idx", "id", "date", "period", "amount", "ln_price", "region", "region_nom", "mrc", "mrc_nom", "code_mun", | |
| 184 | + "mun_nom", "unite_voisinage", "floor_area", "land_area", "age", "lien", "etages", "role_vente", "month", "lat", "lng"]].reset_index(drop=True) | |
| 185 | + T["n_cand"] = [x["n_cand"] for x in infos]; T["asof"] = [x["asof"] for x in infos] | |
| 186 | + T["n_2km_12m"] = nav[:, DS.index(2) * len(TS) + TS.index(12)] | |
| 187 | + T["n_5km_12m"] = nav[:, DS.index(5) * len(TS) + TS.index(12)] | |
| 188 | + save_parquet(T, "targets.parquet") | |
| 189 | + pd.DataFrame({"DS": [DS], "TS": [TS], "KS": [KS]}).to_json(OUT / "grid_meta.json") | |
| 190 | + log(f"rétro-évaluation terminée en {time.time() - t0:.0f} s") | |
added
pipeline/05_analyze.py
+358 −0
@@ -0,0 +1,358 @@ | ||
| 1 | +# -*- coding: utf-8 -*- | |
| 2 | +"""05 — Analyse : métriques par configuration (validation 2024 → choix ; test 2025-26 → verdict), | |
| 3 | +surface Erreur(k, d, t), effets des facteurs (sélection, ajustement, pondération, réconciliation), | |
| 4 | +compétition de modèles M1-M13, hypothèses H1-H5, strates (région, densité, prix), intervalles, | |
| 5 | +règles conditionnelles. Sorties : tables/*.csv, results.json, fig/*.pdf""" | |
| 6 | +from scipy.optimize import nnls | |
| 7 | +import matplotlib; matplotlib.use("Agg") | |
| 8 | +import matplotlib.pyplot as plt | |
| 9 | +from matplotlib import rcParams | |
| 10 | +from common import * | |
| 11 | + | |
| 12 | +# ---------------------------------------------------------------- style (palette de référence, marques fines) | |
| 13 | +PAL = ["#2a78d6", "#eb6834", "#1baf7a", "#eda100", "#e87ba4", "#008300", "#6250d6", "#e34948"] | |
| 14 | +SEQ = ["#cde2fb", "#9ec5f4", "#6da7ec", "#3987e5", "#256abf", "#184f95", "#0d366b"] | |
| 15 | +rcParams.update({"font.family": "Helvetica Neue", "font.size": 9, "axes.edgecolor": "#c3c2b7", "axes.linewidth": 0.6, | |
| 16 | + "axes.grid": True, "grid.color": "#e8e7e3", "grid.linewidth": 0.5, "axes.spines.top": False, | |
| 17 | + "axes.spines.right": False, "legend.frameon": False, "figure.dpi": 150, "savefig.bbox": "tight", | |
| 18 | + "axes.titleweight": "bold", "axes.titlesize": 10, "axes.labelcolor": "#52514e", "xtick.color": "#52514e", | |
| 19 | + "ytick.color": "#52514e", "figure.facecolor": "#fcfcfb", "axes.facecolor": "#fcfcfb"}) | |
| 20 | +from matplotlib.colors import LinearSegmentedColormap | |
| 21 | +CMAP = LinearSegmentedColormap.from_list("seqblue", SEQ) | |
| 22 | +TAB = OUT / "tables"; TAB.mkdir(exist_ok=True); FIG = OUT / "fig" | |
| 23 | +def savefig(fig, name): fig.savefig(FIG / f"{name}.pdf"); fig.savefig(FIG / f"{name}.png"); plt.close(fig) | |
| 24 | + | |
| 25 | +# ---------------------------------------------------------------- chargement | |
| 26 | +T = load_parquet("targets.parquet"); S = load_parquet("sales.parquet").set_index("idx") | |
| 27 | +E1 = pd.read_csv(OUT / "configs_E1.csv"); E2 = pd.read_csv(OUT / "configs_E2.csv") | |
| 28 | +est1 = np.load(OUT / "est_E1.npy"); est2 = np.load(OUT / "est_E2.npy"); sd2 = np.load(OUT / "sd_E2.npy"); nav = np.load(OUT / "n_avail.npy") | |
| 29 | +meta = load_json("index_meta.json"); regions = meta["regions"]; M_EVAL0 = meta["M_EVAL0"]; IDX = np.load(OUT / "index_asof.npy") | |
| 30 | +rid = {r: i for i, r in enumerate(regions)} | |
| 31 | +lnP = T["ln_price"].to_numpy(np.float32) | |
| 32 | +VAL = (T["period"] == "val").to_numpy(); TEST = (T["period"] == "test").to_numpy() | |
| 33 | +log(f"{len(T):,} cibles (val {VAL.sum():,}, test {TEST.sum():,}) ; E1 {est1.shape[1]} ; E2 {est2.shape[1]}") | |
| 34 | + | |
| 35 | +def err_of(est): return (np.exp(est - lnP[:, None]) - 1).astype(np.float32) | |
| 36 | +def metrics(err, sel): | |
| 37 | + e = err[sel]; a = np.abs(e) | |
| 38 | + ok = ~np.isnan(a); n = ok.sum(axis=0) | |
| 39 | + with np.errstate(all="ignore"): | |
| 40 | + return pd.DataFrame({"mdape": np.nanmedian(a, axis=0) * 100, "mape": np.nanmean(a, axis=0) * 100, | |
| 41 | + "biais": np.nanmedian(e, axis=0) * 100, "biais_moy": np.nanmean(e, axis=0) * 100, | |
| 42 | + "p5": np.nanmean(a <= 0.05, axis=0) * 100, "p10": np.nanmean(a <= 0.10, axis=0) * 100, | |
| 43 | + "p20": np.nanmean(a <= 0.20, axis=0) * 100, "couverture": n / sel.sum() * 100, "n": n}) | |
| 44 | +err1 = err_of(est1); err2 = err_of(est2) | |
| 45 | +del est1 | |
| 46 | +res = {} | |
| 47 | + | |
| 48 | +# ---------------------------------------------------------------- densité & strates | |
| 49 | +T["dens"] = pd.cut(T["n_2km_12m"], [-1, 4, 19, 49, 199, 10 ** 9], labels=["< 5", "5-19", "20-49", "50-199", "≥ 200"]) | |
| 50 | +T["prix_q"] = pd.qcut(T["amount"], 5, labels=["Q1 (bas)", "Q2", "Q3", "Q4", "Q5 (haut)"]) | |
| 51 | +T["trimestre"] = pd.to_datetime(T["date"]).dt.to_period("Q").astype(str) | |
| 52 | +T["age_cl"] = pd.cut(T["age"], [-1, 3, 15, 40, 70, 1000], labels=["≤ 3 ans", "4-15", "16-40", "41-70", "> 70"]) | |
| 53 | +T["cluster"] = pd.factorize(T["code_mun"].astype(str) + "|" + T["trimestre"])[0] | |
| 54 | + | |
| 55 | +# ---------------------------------------------------------------- E1 : surface Erreur(k, d, t) | |
| 56 | +mv = metrics(err1, VAL).add_suffix("_val"); mt = metrics(err1, TEST).add_suffix("_test") | |
| 57 | +surf = pd.concat([E1, mv, mt], axis=1) | |
| 58 | +surf.to_csv(TAB / "surface_E1.csv", index=False) | |
| 59 | +ok = surf["couverture_val"] >= 97 | |
| 60 | +best_e1 = surf[ok].sort_values("mdape_val").iloc[0] | |
| 61 | +res["E1_best_val"] = best_e1.to_dict() | |
| 62 | +log(f"E1 meilleure (val, couverture ≥ 97 %) : {best_e1[['sel','d','T','k','adj']].to_dict()} → val {best_e1.mdape_val:.2f} %, test {best_e1.mdape_test:.2f} %") | |
| 63 | +for (s, a), g in surf[ok].groupby(["sel", "adj"]): | |
| 64 | + b = g.sort_values("mdape_val").iloc[0] | |
| 65 | + res.setdefault("E1_best_by_sel_adj", {})[f"{s}|{a}"] = dict(k=int(b.k), d=float(b.d), T=int(b["T"]), mdape_val=float(b.mdape_val), mdape_test=float(b.mdape_test), couverture_test=float(b.couverture_test)) | |
| 66 | + | |
| 67 | +# figure : MdAPE vs k selon d (T = 12, geo, hedo, test) | |
| 68 | +fig, axes = plt.subplots(1, 2, figsize=(9, 3.4), sharey=True) | |
| 69 | +for ax, a in zip(axes, ["none", "hedo"]): | |
| 70 | + for j, d in enumerate([0.5, 1, 2, 5, 10, 25]): | |
| 71 | + g = surf[(surf.sel == "geo") & (surf["T"] == 12) & (surf.adj == a) & (surf.d == d)].sort_values("k") | |
| 72 | + ax.plot(g.k, g.mdape_test, color=PAL[j], lw=2, marker="o", ms=4, label=f"d ≤ {d} km") | |
| 73 | + ax.set_xscale("log"); ax.set_xticks([1, 3, 5, 10, 20, 50]); ax.set_xticklabels([1, 3, 5, 10, 20, 50]) | |
| 74 | + ax.set_xlabel("nombre de comparables k"); ax.set_title({"none": "sans ajustement (temps seulement)", "hedo": "grille hédonique"}[a]) | |
| 75 | +axes[0].set_ylabel("MdAPE, test 2025-26 (%)"); axes[1].legend(ncol=2, fontsize=8) | |
| 76 | +savefig(fig, "f1_mdape_vs_k") | |
| 77 | +# figure : carte d × T au meilleur k (geo, hedo) | |
| 78 | +fig, axes = plt.subplots(1, 2, figsize=(8, 3.2)) | |
| 79 | +for ax, a in zip(axes, ["none", "hedo"]): | |
| 80 | + g = surf[(surf.sel == "geo") & (surf.adj == a)].groupby(["d", "T"])["mdape_test"].min().unstack("T") | |
| 81 | + im = ax.imshow(g.to_numpy(), cmap=CMAP, aspect="auto") | |
| 82 | + ax.set_xticks(range(len(g.columns))); ax.set_xticklabels([f"{t} mois" for t in g.columns]); ax.set_yticks(range(len(g.index))); ax.set_yticklabels([f"{d} km" for d in g.index]) | |
| 83 | + for i in range(g.shape[0]): | |
| 84 | + for j in range(g.shape[1]): | |
| 85 | + v = g.iloc[i, j]; ax.text(j, i, f"{v:.1f}", ha="center", va="center", fontsize=7.5, color="white" if v < np.nanpercentile(g.to_numpy(), 40) else "#0b0b0b") | |
| 86 | + ax.grid(False); ax.set_title({"none": "sans ajustement", "hedo": "grille hédonique"}[a] + " — MdAPE au meilleur k") | |
| 87 | +savefig(fig, "f2_surface_d_T") | |
| 88 | + | |
| 89 | +# ---------------------------------------------------------------- E2 : facteurs | |
| 90 | +mv2 = metrics(err2, VAL).add_suffix("_val"); mt2 = metrics(err2, TEST).add_suffix("_test") | |
| 91 | +fac = pd.concat([E2, mv2, mt2], axis=1); fac.to_csv(TAB / "facteurs_E2.csv", index=False) | |
| 92 | +anchor = (fac.k == 5) & (fac.d == 5) & (fac["T"] == 12) | |
| 93 | +base = anchor & (fac.sel == "geo") & (fac.w == "equal") & (fac.recon == "median") | |
| 94 | +res["E2_adjust_anchor"] = fac[base].set_index("adj")[["mdape_val", "mdape_test", "mape_test", "biais_test", "p10_test", "couverture_test"]].to_dict("index") | |
| 95 | +res["E2_sel_anchor"] = fac[anchor & (fac.adj == "hedo") & (fac.w == "equal") & (fac.recon == "median")].set_index("sel")[["mdape_val", "mdape_test", "p10_test"]].to_dict("index") | |
| 96 | +res["E2_wr_anchor"] = {f"{r.w}|{r.recon}": dict(mdape_val=r.mdape_val, mdape_test=r.mdape_test, p10_test=r.p10_test) | |
| 97 | + for r in fac[anchor & (fac.adj == "hedo") & (fac.sel == "geo")].itertuples()} | |
| 98 | +# effets marginaux moyens (moyenne de MdAPE test sur toutes les autres modalités) | |
| 99 | +marg = {} | |
| 100 | +for f in ["sel", "adj", "w", "recon"]: | |
| 101 | + marg[f] = fac.groupby(f)["mdape_test"].mean().round(2).to_dict() | |
| 102 | +res["E2_marginal"] = marg | |
| 103 | +anc_tab = fac.groupby(["k", "d", "T"])["mdape_test"].min().round(2).to_dict() | |
| 104 | +res["E2_anchor_best"] = {f"{k}|{d}|{t}": v for (k, d, t), v in anc_tab.items()} | |
| 105 | +# figure : ajustements à l'ancre (val et test) | |
| 106 | +fig, ax = plt.subplots(figsize=(7, 3.2)) | |
| 107 | +g = fac[base].set_index("adj").loc[["raw", "none", "conv", "hedo", "hedo_reg", "local", "gbm", "role"]] | |
| 108 | +lab = {"raw": "aucun", "none": "temps seul", "conv": "grille conv.", "hedo": "hédonique", "hedo_reg": "hédo. régionale", "local": "locale", "gbm": "GBM", "role": "ratio au rôle"} | |
| 109 | +x = np.arange(len(g)); ax.bar(x - 0.2, g.mdape_val, 0.38, color=PAL[0], label="validation 2024"); ax.bar(x + 0.2, g.mdape_test, 0.38, color=PAL[1], label="test 2025-26") | |
| 110 | +for i, (v1, v2) in enumerate(zip(g.mdape_val, g.mdape_test)): | |
| 111 | + ax.text(i - 0.2, v1 + 0.2, f"{v1:.1f}", ha="center", fontsize=7); ax.text(i + 0.2, v2 + 0.2, f"{v2:.1f}", ha="center", fontsize=7) | |
| 112 | +ax.set_xticks(x); ax.set_xticklabels([lab[a] for a in g.index]); ax.set_ylabel("MdAPE (%)"); ax.legend(); ax.set_title("Type d'ajustement — 5 comparables les plus proches, ≤ 5 km, ≤ 12 mois, médiane") | |
| 113 | +savefig(fig, "f3_ajustements") | |
| 114 | + | |
| 115 | +# ---------------------------------------------------------------- modèles de référence (M10, M11) et compétition | |
| 116 | +def level_shift(m_asof): | |
| 117 | + mt_ = T["month"].to_numpy(); r = T["region"].map(rid).to_numpy() | |
| 118 | + return IDX[mt_ - M_EVAL0, r, mt_ - 1] - IDX[mt_ - M_EVAL0, r, m_asof - 1] | |
| 119 | +Sx = S.loc[T["idx"]] | |
| 120 | +XF = np.column_stack([Sx.ln_floor, Sx.ln_land, Sx.age, Sx.age ** 2 / 100, Sx.lien == 1, Sx.lien >= 2, Sx.genre == 1, Sx.genre == 2, Sx.genre == 3, Sx.etages, Sx.neuf]).astype(float) | |
| 121 | +M10 = np.full(len(T), np.nan); M11 = np.full(len(T), np.nan); M11b = np.full(len(T), np.nan) | |
| 122 | +for a in load_json("asof.json")["asof"]: | |
| 123 | + mask = (T["asof"] == a).to_numpy(); period = str(a) | |
| 124 | + if not mask.any(): continue | |
| 125 | + mod = load_json(f"models_{period}.json"); beta = np.array([mod["beta_global"][f] for f in mod["feats"]]) | |
| 126 | + shift = level_shift(mod["m_asof"]) | |
| 127 | + fe = (T["region"].map(mod["region_fe"]).fillna(0) + T["mrc"].map(mod["mrc_eff"]).fillna(0) + T["code_mun"].map(mod["mun_eff"]).fillna(0) | |
| 128 | + + (T["code_mun"].astype(str) + "|" + T["unite_voisinage"].astype(str)).map(mod["uv_eff"]).fillna(0)).to_numpy() | |
| 129 | + M10[mask] = (XF @ beta + fe + shift)[mask] | |
| 130 | + M11[mask] = (np.load(OUT / f"gbm_role_pred_{period}.npy")[T["idx"].to_numpy()] + shift)[mask] | |
| 131 | + M11b[mask] = (np.load(OUT / f"gbm_pred_{period}.npy")[T["idx"].to_numpy()] + shift)[mask] | |
| 132 | + | |
| 133 | +def e2col(sel, k, d, Tm, adj, w, recon): | |
| 134 | + m = (E2.sel == sel) & (E2.k == k) & (E2.d == d) & (E2["T"] == Tm) & (E2.adj == adj) & (E2.w == w) & (E2.recon == recon) | |
| 135 | + return int(E2.index[m][0]) | |
| 136 | +def e1col(sel, k, d, Tm, adj): | |
| 137 | + m = (E1.sel == sel) & (E1.k == k) & (E1.d == d) & (E1["T"] == Tm) & (E1.adj == adj) | |
| 138 | + return int(E1.index[m][0]) | |
| 139 | +MODELS = { | |
| 140 | + "M1 médiane des 3 plus proches (temps seul)": ("E2", e2col("geo", 3, 2, 6, "none", "equal", "median")), | |
| 141 | + "M2 moyenne des 5 plus proches (temps seul)": ("E2", e2col("geo", 5, 5, 12, "none", "equal", "mean")), | |
| 142 | + "M3 sélection par similarité physique": ("E2", e2col("sim", 5, 5, 12, "none", "equal", "median")), | |
| 143 | + "M4 pondération par distance": ("E2", e2col("geo", 5, 5, 12, "none", "dist", "mean")), | |
| 144 | + "M5 pondération par récence": ("E2", e2col("geo", 5, 5, 12, "none", "rec", "mean")), | |
| 145 | + "M6 grille d'ajustement conventionnelle": ("E2", e2col("geo", 5, 5, 12, "conv", "equal", "median")), | |
| 146 | + "M7 grille d'ajustement hédonique": ("E2", e2col("geo", 5, 5, 12, "hedo", "equal", "median")), | |
| 147 | + "M8a grille locale (régression sur le voisinage)": ("E2", e2col("geo", 5, 5, 12, "local", "equal", "median")), | |
| 148 | + "M8b ajustements non linéaires (GBM)": ("E2", e2col("geo", 5, 5, 12, "gbm", "equal", "median")), | |
| 149 | + "M13 ratio au rôle (5 plus proches, même rôle)": ("E2", e2col("geo", 5, 5, 12, "role", "equal", "median")), | |
| 150 | + "M10 hédonique globale directe": ("X", M10), | |
| 151 | + "M11a gradient boosting (sans rôle)": ("X", M11b), | |
| 152 | + "M11b gradient boosting (avec rôle)": ("X", M11), | |
| 153 | +} | |
| 154 | +def est_of(spec): | |
| 155 | + kind, c = spec | |
| 156 | + return est2[:, c].astype(float) if kind == "E2" else c | |
| 157 | +# M0 : meilleure configuration empirique choisie sur la validation (E1 ∪ E2, couverture ≥ 90 %) | |
| 158 | +cand_all = pd.concat([surf.assign(src="E1"), fac.assign(src="E2")], ignore_index=True) | |
| 159 | +cand_ok = cand_all[cand_all.couverture_val >= 97].sort_values("mdape_val") | |
| 160 | +b0 = cand_ok.iloc[0] | |
| 161 | +M0_est = (est2[:, int(b0.cfg)] if b0.src == "E2" else np.load(OUT / "est_E1.npy", mmap_mode="r")[:, int(b0.cfg)]).astype(float) | |
| 162 | +res["M0"] = b0[["src", "sel", "d", "T", "k", "adj", "w", "recon", "mdape_val", "mdape_test", "couverture_test"]].to_dict() | |
| 163 | +res["M0_sans_contrainte"] = cand_all[cand_all.couverture_val >= 50].sort_values("mdape_val").iloc[0][["src", "sel", "d", "T", "k", "adj", "w", "recon", "mdape_val", "mdape_test", "couverture_val", "couverture_test"]].to_dict() | |
| 164 | +MODELS["M0 meilleure configuration empirique (choisie en 2024, couverture ≥ 97 %)"] = ("X", M0_est) | |
| 165 | +# M9 : réconciliation optimisée = combinaison convexe (NNLS sur la validation) d'estimations de comparables | |
| 166 | +stack_cols = {a: e2col("geo", 5, 5, 12, a, "equal", "median") for a in ["none", "conv", "hedo", "hedo_reg", "local", "gbm", "role"]} | |
| 167 | +def stack(names, mats, fit_mask): | |
| 168 | + """Combinaison convexe apprise par NNLS sur la validation ; si une composante manque pour une cible, | |
| 169 | + les poids des composantes disponibles sont renormalisés.""" | |
| 170 | + A = np.column_stack(mats); okr = ~np.isnan(A).any(axis=1) | |
| 171 | + w, _ = nnls(A[okr & fit_mask], lnP[okr & fit_mask].astype(float)); w = w / w.sum() | |
| 172 | + W = np.where(np.isnan(A), 0.0, w[None, :]); den = W.sum(axis=1) | |
| 173 | + out = np.nansum(A * W, axis=1) / np.where(den > 0, den, np.nan) | |
| 174 | + return out, dict(zip(names, w.round(3))) | |
| 175 | +M9_est, w9 = stack(list(stack_cols), [est2[:, c].astype(float) for c in stack_cols.values()], VAL) | |
| 176 | +res["M9_weights"] = w9; MODELS["M9 réconciliation optimisée (poids appris en 2024)"] = ("X", M9_est) | |
| 177 | +M12_est, w12 = stack(["gbm_role", "comps_gbm", "comps_hedo"], [M11, est2[:, stack_cols["gbm"]].astype(float), est2[:, stack_cols["hedo"]].astype(float)], VAL) | |
| 178 | +res["M12_weights"] = w12; MODELS["M12 hybride comparables + GBM (poids appris en 2024)"] = ("X", M12_est) | |
| 179 | + | |
| 180 | +rows = [] | |
| 181 | +ERR = {} | |
| 182 | +for name, spec in MODELS.items(): | |
| 183 | + e = (np.exp(est_of(spec) - lnP) - 1).astype(np.float32); ERR[name] = e | |
| 184 | + for per, mask in [("val", VAL), ("test", TEST)]: | |
| 185 | + m = metrics(e[:, None], mask).iloc[0].to_dict(); m.update(modele=name, periode=per); rows.append(m) | |
| 186 | +comp = pd.DataFrame(rows); comp.to_csv(TAB / "competition_modeles.csv", index=False) | |
| 187 | +pd.DataFrame({n.split(" ")[0]: est_of(s_) for n, s_ in MODELS.items()}).assign(idx=T["idx"].to_numpy()).to_parquet(OUT / "models_est.parquet", index=False) | |
| 188 | +res["competition"] = comp[comp.periode == "test"].set_index("modele")[["mdape", "mape", "biais", "p5", "p10", "p20", "couverture"]].round(2).to_dict("index") | |
| 189 | +res["competition_val"] = comp[comp.periode == "val"].set_index("modele")[["mdape", "p10", "couverture"]].round(2).to_dict("index") | |
| 190 | +log("compétition (test) :\n" + comp[comp.periode == "test"][["modele", "mdape", "mape", "biais", "p10", "couverture"]].round(2).to_string(index=False)) | |
| 191 | +fig, ax = plt.subplots(figsize=(7.5, 5)) | |
| 192 | +g = comp[comp.periode == "test"].sort_values("mdape", ascending=False) | |
| 193 | +cols = [PAL[0] if ("M1" in n and "M1 " not in n) or n.startswith(("M10", "M11", "M12")) else PAL[2] if n.startswith(("M0", "M9")) else PAL[1] for n in g.modele] | |
| 194 | +cols = [PAL[0] if n.startswith(("M10", "M11", "M12")) else (PAL[2] if n.startswith(("M0", "M9")) else PAL[1]) for n in g.modele] | |
| 195 | +ax.barh(g.modele, g.mdape, color=cols, height=0.65) | |
| 196 | +for i, v in enumerate(g.mdape): ax.text(v + 0.15, i, f"{v:.1f} %", va="center", fontsize=8) | |
| 197 | +ax.set_xlabel("MdAPE, test 2025-26 (%)"); ax.set_title("Compétition de modèles (orange : comparables ; vert : optimisés ; bleu : références externes)") | |
| 198 | +savefig(fig, "f4_competition") | |
| 199 | + | |
| 200 | +# ---------------------------------------------------------------- bootstrap par grappes (municipalité × trimestre) des différences appariées | |
| 201 | +rng = np.random.default_rng(7) | |
| 202 | +def paired(eA, eB, mask, B=400): | |
| 203 | + a = np.abs(eA[mask]); b = np.abs(eB[mask]); cl = T["cluster"].to_numpy()[mask] | |
| 204 | + okp = ~np.isnan(a) & ~np.isnan(b); a, b, cl = a[okp], b[okp], cl[okp] | |
| 205 | + d = a - b | |
| 206 | + cl_ids, inv = np.unique(cl, return_inverse=True); nc = len(cl_ids) | |
| 207 | + sums = np.bincount(inv, d, nc); cnts = np.bincount(inv, minlength=nc) | |
| 208 | + boots_mean = []; boots_med = [] | |
| 209 | + order = np.argsort(inv); d_sorted = d[order]; starts = np.searchsorted(inv[order], np.arange(nc)); ends = np.append(starts[1:], len(d)) | |
| 210 | + for _ in range(B): | |
| 211 | + pick = rng.integers(0, nc, nc) | |
| 212 | + boots_mean.append(sums[pick].sum() / cnts[pick].sum()) | |
| 213 | + if _ < 150: | |
| 214 | + idxs = np.concatenate([np.arange(starts[p], ends[p]) for p in pick]) | |
| 215 | + boots_med.append(np.median(a[order][idxs]) - np.median(b[order][idxs])) | |
| 216 | + return dict(n=int(len(d)), diff_mape=float(d.mean() * 100), ci_mape=[float(np.percentile(boots_mean, 2.5) * 100), float(np.percentile(boots_mean, 97.5) * 100)], | |
| 217 | + diff_mdape=float((np.median(a) - np.median(b)) * 100), ci_mdape=[float(np.percentile(boots_med, 2.5) * 100), float(np.percentile(boots_med, 97.5) * 100)], | |
| 218 | + part_A_meilleur=float((a < b).mean() * 100)) | |
| 219 | +PAIRS = [ | |
| 220 | + ("M7 grille d'ajustement hédonique", "M2 moyenne des 5 plus proches (temps seul)"), | |
| 221 | + ("M7 grille d'ajustement hédonique", "M6 grille d'ajustement conventionnelle"), | |
| 222 | + ("M8a grille locale (régression sur le voisinage)", "M7 grille d'ajustement hédonique"), | |
| 223 | + ("M8b ajustements non linéaires (GBM)", "M7 grille d'ajustement hédonique"), | |
| 224 | + ("M11b gradient boosting (avec rôle)", "M0 meilleure configuration empirique (choisie en 2024, couverture ≥ 97 %)"), | |
| 225 | + ("M11a gradient boosting (sans rôle)", "M0 meilleure configuration empirique (choisie en 2024, couverture ≥ 97 %)"), | |
| 226 | + ("M12 hybride comparables + GBM (poids appris en 2024)", "M11b gradient boosting (avec rôle)"), | |
| 227 | + ("M13 ratio au rôle (5 plus proches, même rôle)", "M7 grille d'ajustement hédonique"), | |
| 228 | + ("M3 sélection par similarité physique", "M2 moyenne des 5 plus proches (temps seul)"), | |
| 229 | + ("M9 réconciliation optimisée (poids appris en 2024)", "M7 grille d'ajustement hédonique"), | |
| 230 | + ("M1 médiane des 3 plus proches (temps seul)", "M2 moyenne des 5 plus proches (temps seul)"), | |
| 231 | +] | |
| 232 | +res["paires"] = {f"{a} vs {b}": paired(ERR[a], ERR[b], TEST) for a, b in PAIRS} | |
| 233 | +pd.DataFrame([dict(A=a, B=b, **res["paires"][f"{a} vs {b}"]) for a, b in PAIRS]).to_csv(TAB / "paires_bootstrap.csv", index=False) | |
| 234 | + | |
| 235 | +# ---------------------------------------------------------------- strates | |
| 236 | +def strata(err, by): | |
| 237 | + out = [] | |
| 238 | + for lev, g in T[TEST].groupby(by, observed=True): | |
| 239 | + m = np.zeros(len(T), bool); m[g.index] = True | |
| 240 | + r = metrics(err[:, None], m).iloc[0].to_dict(); r[by] = str(lev); out.append(r) | |
| 241 | + return pd.DataFrame(out) | |
| 242 | +KEY = ["M2 moyenne des 5 plus proches (temps seul)", "M7 grille d'ajustement hédonique", "M8a grille locale (régression sur le voisinage)", | |
| 243 | + "M13 ratio au rôle (5 plus proches, même rôle)", "M0 meilleure configuration empirique (choisie en 2024, couverture ≥ 97 %)", "M10 hédonique globale directe", "M11b gradient boosting (avec rôle)", "M12 hybride comparables + GBM (poids appris en 2024)"] | |
| 244 | +for by in ["region_nom", "dens", "prix_q", "trimestre", "age_cl", "lien"]: | |
| 245 | + tabs = [] | |
| 246 | + for name in KEY: | |
| 247 | + s = strata(ERR[name], by); s["modele"] = name.split(" ")[0]; tabs.append(s) | |
| 248 | + st = pd.concat(tabs); st.to_csv(TAB / f"strates_{by}.csv", index=False) | |
| 249 | + res[f"strates_{by}"] = st.pivot(index=by, columns="modele", values="mdape").round(2).to_dict("index") | |
| 250 | + res[f"strates_{by}_n"] = st[st.modele == "M7"].set_index(by)["n"].to_dict() | |
| 251 | +# figure : densité | |
| 252 | +fig, ax = plt.subplots(figsize=(6.5, 3.4)) | |
| 253 | +st = pd.read_csv(TAB / "strates_dens.csv") | |
| 254 | +order = ["< 5", "5-19", "20-49", "50-199", "≥ 200"] | |
| 255 | +for j, (name, lab_) in enumerate([("M2", "M2 5 plus proches"), ("M7", "M7 grille hédonique"), ("M0", "M0 meilleure config."), ("M11b", "M11b GBM"), ("M12", "M12 hybride")]): | |
| 256 | + g = st[st.modele == name].set_index("dens").loc[order] | |
| 257 | + ax.plot(order, g.mdape, color=PAL[j], lw=2, marker="o", ms=4, label=lab_) | |
| 258 | +ax.set_xlabel("ventes disponibles dans un rayon de 2 km au cours des 12 mois précédents"); ax.set_ylabel("MdAPE, test (%)"); ax.legend(fontsize=8); ax.set_title("Précision selon la liquidité du sous-marché") | |
| 259 | +savefig(fig, "f5_densite") | |
| 260 | + | |
| 261 | +# ---------------------------------------------------------------- hypothèses | |
| 262 | +H = {} | |
| 263 | +# H1 : proche mais dissemblable vs loin mais semblable (k = 5, T = 12, sans ajustement) | |
| 264 | +H["H1"] = {"geo_d1": surf.loc[e1col("geo", 5, 1, 12, "none"), ["mdape_test", "couverture_test"]].to_dict(), | |
| 265 | + "sim_d10": surf.loc[e1col("sim", 5, 10, 12, "none"), ["mdape_test", "couverture_test"]].to_dict(), | |
| 266 | + "geo_d10": surf.loc[e1col("geo", 5, 10, 12, "none"), ["mdape_test", "couverture_test"]].to_dict(), | |
| 267 | + "mix_d10": surf.loc[e1col("mix", 5, 10, 12, "none"), ["mdape_test", "couverture_test"]].to_dict(), | |
| 268 | + "sim_d10_hedo": surf.loc[e1col("sim", 5, 10, 12, "hedo"), ["mdape_test", "couverture_test"]].to_dict(), | |
| 269 | + "geo_d10_hedo": surf.loc[e1col("geo", 5, 10, 12, "hedo"), ["mdape_test", "couverture_test"]].to_dict(), | |
| 270 | + "mix_d10_hedo": surf.loc[e1col("mix", 5, 10, 12, "hedo"), ["mdape_test", "couverture_test"]].to_dict()} | |
| 271 | +# H2 : k optimal selon la densité (geo, d = 10, T = 12, hedo et none) | |
| 272 | +h2 = {} | |
| 273 | +for a in ["none", "hedo"]: | |
| 274 | + for lev in order: | |
| 275 | + m = TEST & (T["dens"] == lev).to_numpy() | |
| 276 | + cols_ = [e1col("geo", k, 10, 12, a) for k in [1, 3, 5, 7, 10, 15, 20, 30, 50]] | |
| 277 | + md = np.nanmedian(np.abs(err1[m][:, cols_]), axis=0) * 100 | |
| 278 | + h2[f"{a}|{lev}"] = dict(zip([1, 3, 5, 7, 10, 15, 20, 30, 50], md.round(2))) | |
| 279 | +H["H2"] = h2 | |
| 280 | +fig, axes = plt.subplots(1, 2, figsize=(9, 3.3), sharey=True) | |
| 281 | +for ax, a in zip(axes, ["none", "hedo"]): | |
| 282 | + for j, lev in enumerate(order): | |
| 283 | + v = h2[f"{a}|{lev}"]; ax.plot(list(v), list(v.values()), color=PAL[j], lw=2, marker="o", ms=3.5, label=f"densité {lev}") | |
| 284 | + ax.set_xscale("log"); ax.set_xticks([1, 3, 5, 10, 20, 50]); ax.set_xticklabels([1, 3, 5, 10, 20, 50]); ax.set_xlabel("k"); ax.set_title({"none": "temps seul", "hedo": "grille hédonique"}[a] + " — d ≤ 10 km, 12 mois") | |
| 285 | +axes[0].set_ylabel("MdAPE, test (%)"); axes[1].legend(fontsize=8) | |
| 286 | +savefig(fig, "f6_k_par_densite") | |
| 287 | +# H3 : local vs hédonique fixe selon la densité, aux trois ancres | |
| 288 | +h3 = {} | |
| 289 | +for (k, d, Tm) in [(3, 2, 6), (5, 5, 12), (10, 10, 24)]: | |
| 290 | + for lev in order: | |
| 291 | + m = TEST & (T["dens"] == lev).to_numpy() | |
| 292 | + h3[f"{k}|{d}|{Tm}|{lev}"] = {a: float(np.nanmedian(np.abs(err2[m][:, e2col("geo", k, d, Tm, a, "equal", "median")])) * 100) for a in ["conv", "hedo", "hedo_reg", "local", "gbm"]} | |
| 293 | +H["H3"] = h3 | |
| 294 | +# H4 : pondération selon la densité (ancre 5/5/12, hedo) | |
| 295 | +h4 = {} | |
| 296 | +for lev in order: | |
| 297 | + m = TEST & (T["dens"] == lev).to_numpy() | |
| 298 | + h4[lev] = {f"{w}|{r}": float(np.nanmedian(np.abs(err2[m][:, e2col("geo", 5, 5, 12, "hedo", w, r)])) * 100) for (w, r) in [("equal", "median"), ("equal", "mean"), ("dist", "mean"), ("rec", "mean"), ("sim", "mean"), ("all", "mean"), ("all", "median")]} | |
| 299 | +H["H4"] = h4 | |
| 300 | +# H5 : méthodes simples vs ML selon l'homogénéité (dispersion des prix ajustés des 10 voisins à 10 km / 24 mois) | |
| 301 | +het = sd2[:, e2col("geo", 10, 10, 24, "hedo", "equal", "median")].astype(float) | |
| 302 | +T["heterog"] = pd.qcut(pd.Series(np.where(np.isnan(het), np.nanmax(het), het)), 3, labels=["homogène", "moyen", "hétérogène"]) | |
| 303 | +h5 = {} | |
| 304 | +for lev in ["homogène", "moyen", "hétérogène"]: | |
| 305 | + m = TEST & (T["heterog"] == lev).to_numpy() | |
| 306 | + h5[lev] = {n.split(" ")[0]: float(np.nanmedian(np.abs(ERR[n][m])) * 100) for n in KEY} | |
| 307 | + h5[lev]["sd_median"] = float(np.nanmedian(het[m])) | |
| 308 | +H["H5"] = h5 | |
| 309 | +res["hypotheses"] = H | |
| 310 | + | |
| 311 | +# ---------------------------------------------------------------- intervalles d'incertitude (couverture) | |
| 312 | +intervals = {} | |
| 313 | +for a in ["none", "hedo", "gbm"]: | |
| 314 | + for (k, d, Tm) in [(5, 5, 12), (10, 10, 24)]: | |
| 315 | + c = e2col("geo", k, d, Tm, a, "equal", "median"); e = est2[:, c].astype(float); s = sd2[:, c].astype(float) | |
| 316 | + okv = VAL & ~np.isnan(e) & ~np.isnan(s) & (s > 0) | |
| 317 | + z_needed = np.abs(lnP[okv] - e[okv]) / s[okv]; z80 = float(np.percentile(z_needed, 80)) | |
| 318 | + okt = TEST & ~np.isnan(e) & ~np.isnan(s) & (s > 0) | |
| 319 | + cov_nominal = float(np.mean(np.abs(lnP[okt] - e[okt]) <= 1.2816 * s[okt]) * 100) | |
| 320 | + cov_cal = float(np.mean(np.abs(lnP[okt] - e[okt]) <= z80 * s[okt]) * 100) | |
| 321 | + width = float(np.median(2 * z80 * s[okt]) * 100) | |
| 322 | + intervals[f"{a}|{k}|{d}|{Tm}"] = dict(z80_val=z80, couverture_nominale_80=cov_nominal, couverture_calibree=cov_cal, largeur_mediane_pct=width, n_test=int(okt.sum())) | |
| 323 | +res["intervalles"] = intervals | |
| 324 | + | |
| 325 | +# ---------------------------------------------------------------- règles conditionnelles (système de recommandation) | |
| 326 | +cond = {}; est_cond = np.full(len(T), np.nan) | |
| 327 | +for lev in order: | |
| 328 | + mv_ = VAL & (T["dens"] == lev).to_numpy(); mt_ = TEST & (T["dens"] == lev).to_numpy() | |
| 329 | + m1 = metrics(err1, mv_); m2 = metrics(err2, mv_) | |
| 330 | + allm = pd.concat([pd.concat([E1, m1], axis=1).assign(src="E1"), pd.concat([E2, m2], axis=1).assign(src="E2")], ignore_index=True) | |
| 331 | + b = allm[allm.couverture >= 95].sort_values("mdape").iloc[0] | |
| 332 | + e = (err1 if b.src == "E1" else err2)[:, int(b.cfg)] | |
| 333 | + est_cond[mt_] = e[mt_] | |
| 334 | + cond[lev] = dict(b[["src", "sel", "d", "T", "k", "adj", "w", "recon", "mdape", "couverture"]].to_dict(), mdape_test=float(np.nanmedian(np.abs(e[mt_])) * 100), n_test=int(mt_.sum())) | |
| 335 | +res["regles_conditionnelles"] = cond | |
| 336 | +res["regles_conditionnelles_global_test"] = dict(mdape=float(np.nanmedian(np.abs(est_cond[TEST])) * 100), p10=float(np.nanmean(np.abs(est_cond[TEST]) <= 0.1) * 100), couverture=float(np.mean(~np.isnan(est_cond[TEST])) * 100)) | |
| 337 | +pd.DataFrame(cond).T.to_csv(TAB / "regles_conditionnelles.csv") | |
| 338 | + | |
| 339 | +# ---------------------------------------------------------------- stabilité temporelle, indice, divers | |
| 340 | +tq = [] | |
| 341 | +for q, g in T[TEST].groupby("trimestre"): | |
| 342 | + m = np.zeros(len(T), bool); m[g.index] = True | |
| 343 | + tq.append(dict(trimestre=q, n=int(m.sum()), **{n.split(" ")[0]: float(np.nanmedian(np.abs(ERR[n][m])) * 100) for n in KEY}, | |
| 344 | + **{"biais_" + n.split(" ")[0]: float(np.nanmedian(ERR[n][m]) * 100) for n in KEY})) | |
| 345 | +pd.DataFrame(tq).to_csv(TAB / "stabilite_trimestres.csv", index=False); res["stabilite"] = tq | |
| 346 | +fig, ax = plt.subplots(figsize=(7, 3.2)) | |
| 347 | +idxf = pd.read_csv(OUT / "index_full.csv") | |
| 348 | +for j, (r, lab_) in enumerate([("QC", "Québec"), ("06", "Montréal"), ("03", "Capitale-Nationale"), ("07", "Outaouais"), ("05", "Estrie")]): | |
| 349 | + g = idxf[idxf.region.astype(str).str.zfill(2) == r] if r != "QC" else idxf[idxf.region == "QC"] | |
| 350 | + if len(g): ax.plot(pd.to_datetime(g.ym), np.exp(g.log_idx) * 100, color=PAL[j], lw=2 if r == "QC" else 1.4, label=lab_) | |
| 351 | +ax.set_ylabel("indice hédonique (2021-01 = 100)"); ax.legend(fontsize=8); ax.set_title("Indice de prix mensuel (plein échantillon, descriptif)") | |
| 352 | +savefig(fig, "f7_indice") | |
| 353 | +res["dataset"] = load_json("dataset_stats.json") | |
| 354 | +res["n_targets"] = dict(val=int(VAL.sum()), test=int(TEST.sum())) | |
| 355 | +res["densite_distribution"] = T[TEST]["dens"].value_counts().sort_index().to_dict() | |
| 356 | +res["n_cand_median"] = dict(val=float(np.median(T[VAL].n_cand)), test=float(np.median(T[TEST].n_cand))) | |
| 357 | +save_json(res, "results.json") | |
| 358 | +log("analyse terminée") | |
added
pipeline/06_report.py
+313 −0
@@ -0,0 +1,313 @@ | ||
| 1 | +# -*- coding: utf-8 -*- | |
| 2 | +"""06 — Rapport : assemble docs/RAPPORT.md à partir de results.json et des tables (texte + tableaux + figures). | |
| 3 | +Les phrases d'interprétation sont calculées à partir des résultats (jamais de chiffre écrit à la main).""" | |
| 4 | +from common import * | |
| 5 | + | |
| 6 | +R = load_json("results.json"); TAB = OUT / "tables" | |
| 7 | +DOCS = ROOT / "docs"; DOCS.mkdir(exist_ok=True) | |
| 8 | +D = R["dataset"]; comp = R["competition"]; compv = R["competition_val"] | |
| 9 | +def f1(x): return f"{x:.1f}".replace(".", ",") | |
| 10 | +def f2(x): return f"{x:.2f}".replace(".", ",") | |
| 11 | +def n_(x): return f"{int(round(x)):,}".replace(",", " ") | |
| 12 | +def pct(x): return f"{f1(x)} %" | |
| 13 | +def md_table(df, cols=None, fmt=None, index=False): | |
| 14 | + if cols: df = df[cols] | |
| 15 | + if index: df = df.reset_index() | |
| 16 | + fmt = fmt or {} | |
| 17 | + head = "| " + " | ".join(str(c) for c in df.columns) + " |\n|" + "|".join("---" for _ in df.columns) + "|\n" | |
| 18 | + rows = [] | |
| 19 | + for _, r in df.iterrows(): | |
| 20 | + cells = [] | |
| 21 | + for c in df.columns: | |
| 22 | + v = r[c] | |
| 23 | + if c in fmt: cells.append(fmt[c](v)) | |
| 24 | + elif isinstance(v, float): cells.append(f1(v)) | |
| 25 | + else: cells.append(str(v)) | |
| 26 | + rows.append("| " + " | ".join(cells) + " |") | |
| 27 | + return head + "\n".join(rows) + "\n" | |
| 28 | +short = lambda name: name.split(" ")[0] | |
| 29 | +LAB = {"raw": "aucun ajustement", "none": "temps seulement", "conv": "grille conventionnelle", "hedo": "grille hédonique", "hedo_reg": "grille hédonique régionale", | |
| 30 | + "local": "régression locale", "gbm": "ajustements GBM", "role": "ratio au rôle", "geo": "les plus proches (distance)", "sim": "les plus semblables (physique)", | |
| 31 | + "mix": "compromis distance / similarité", "adj": "les moins ajustés", "equal": "égale", "dist": "1/distance", "rec": "récence", "all": "distance × récence × similarité", | |
| 32 | + "median": "médiane", "mean": "moyenne"} | |
| 33 | + | |
| 34 | +# ------------------------------------------------------------------ tableaux | |
| 35 | +excl = pd.DataFrame(D["exclusions"], columns=["Règle", "n"]); excl["n"] = excl["n"].map(n_) | |
| 36 | +cm = pd.DataFrame(comp).T; cm.index.name = "Modèle"; cm = cm.reset_index() | |
| 37 | +cm["val"] = cm["Modèle"].map(lambda m: compv[m]["mdape"]) | |
| 38 | +cm = cm.sort_values("mdape") | |
| 39 | +t_comp = md_table(cm.rename(columns={"mdape": "MdAPE test", "val": "MdAPE val.", "mape": "MAPE", "biais": "Biais (médian)", "p5": "±5 %", "p10": "±10 %", "p20": "±20 %", "couverture": "Couverture"}), | |
| 40 | + ["Modèle", "MdAPE val.", "MdAPE test", "MAPE", "Biais (médian)", "±5 %", "±10 %", "±20 %", "Couverture"]) | |
| 41 | +surf = pd.read_csv(TAB / "surface_E1.csv") | |
| 42 | +def surf_tab(sel, adj, T_): | |
| 43 | + g = surf[(surf.sel == sel) & (surf.adj == adj) & (surf["T"] == T_)].pivot(index="d", columns="k", values="mdape_test") | |
| 44 | + g.index = [f"≤ {d} km" for d in g.index]; g.columns = [f"k = {k}" for k in g.columns]; g.index.name = "rayon" | |
| 45 | + return md_table(g.round(1), index=True) | |
| 46 | +kb = surf[(surf.sel == "geo") & (surf.adj == "hedo") & (surf.d == 10) & (surf["T"] == 12)].sort_values("k") | |
| 47 | +fac = pd.read_csv(TAB / "facteurs_E2.csv") | |
| 48 | +anchor = fac[(fac.k == 5) & (fac.d == 5) & (fac["T"] == 12)] | |
| 49 | +t_adj = anchor[(anchor.sel == "geo") & (anchor.w == "equal") & (anchor.recon == "median")].set_index("adj").loc[["raw", "none", "conv", "hedo", "hedo_reg", "local", "gbm", "role"]] | |
| 50 | +t_adj = t_adj.reset_index(); t_adj["Ajustement"] = t_adj["adj"].map(LAB) | |
| 51 | +t_adj_md = md_table(t_adj.rename(columns={"mdape_val": "MdAPE val.", "mdape_test": "MdAPE test", "mape_test": "MAPE test", "biais_test": "Biais", "p10_test": "±10 %", "couverture_test": "Couverture"}), | |
| 52 | + ["Ajustement", "MdAPE val.", "MdAPE test", "MAPE test", "Biais", "±10 %", "Couverture"]) | |
| 53 | +t_sel = anchor[(anchor.adj == "hedo") & (anchor.w == "equal") & (anchor.recon == "median")].copy(); t_sel["Sélection"] = t_sel["sel"].map(LAB) | |
| 54 | +t_sel_md = md_table(t_sel.rename(columns={"mdape_val": "MdAPE val.", "mdape_test": "MdAPE test", "p10_test": "±10 %", "biais_test": "Biais"}), ["Sélection", "MdAPE val.", "MdAPE test", "Biais", "±10 %"]) | |
| 55 | +t_wr = anchor[(anchor.adj == "hedo") & (anchor.sel == "geo")].copy(); t_wr["Pondération"] = t_wr["w"].map(LAB); t_wr["Réconciliation"] = t_wr["recon"].map(LAB) | |
| 56 | +t_wr_md = md_table(t_wr.rename(columns={"mdape_val": "MdAPE val.", "mdape_test": "MdAPE test", "p10_test": "±10 %", "biais_test": "Biais"}), ["Pondération", "Réconciliation", "MdAPE val.", "MdAPE test", "Biais", "±10 %"]) | |
| 57 | +pairs = pd.read_csv(TAB / "paires_bootstrap.csv") | |
| 58 | +pairs["Comparaison"] = pairs.apply(lambda r: f"{short(r.A)} − {short(r.B)}", axis=1) | |
| 59 | +pairs["Δ MdAPE (pts) [IC 95 %]"] = pairs.apply(lambda r: f"{f2(r.diff_mdape)} [{f2(eval(r.ci_mdape)[0])} ; {f2(eval(r.ci_mdape)[1])}]", axis=1) | |
| 60 | +pairs["Δ MAPE (pts) [IC 95 %]"] = pairs.apply(lambda r: f"{f2(r.diff_mape)} [{f2(eval(r.ci_mape)[0])} ; {f2(eval(r.ci_mape)[1])}]", axis=1) | |
| 61 | +pairs["A plus précis (% des cibles)"] = pairs["part_A_meilleur"].map(f1) | |
| 62 | +t_pairs = md_table(pairs, ["Comparaison", "Δ MdAPE (pts) [IC 95 %]", "Δ MAPE (pts) [IC 95 %]", "A plus précis (% des cibles)", "n"], fmt={"n": n_}) | |
| 63 | +def strata_md(by, label): | |
| 64 | + st = pd.read_csv(TAB / f"strates_{by}.csv") | |
| 65 | + g = st.pivot(index=by, columns="modele", values="mdape").round(1) | |
| 66 | + g["n"] = st[st.modele == "M7"].set_index(by)["n"] | |
| 67 | + g = g[["n", "M2", "M7", "M8a", "M13", "M0", "M10", "M11b", "M12"]] | |
| 68 | + if by == "dens": g = g.loc[["< 5", "5-19", "20-49", "50-199", "≥ 200"]] | |
| 69 | + if by == "region_nom": g = g.sort_values("n", ascending=False) | |
| 70 | + g.index.name = label | |
| 71 | + return md_table(g, index=True, fmt={"n": n_}) | |
| 72 | +H = R["hypotheses"] | |
| 73 | +h2 = pd.DataFrame({k.split("|")[1]: v for k, v in H["H2"].items() if k.startswith("hedo")}).T; h2.index.name = "densité"; h2.columns = [f"k = {c}" for c in h2.columns] | |
| 74 | +h2_md = md_table(h2.loc[["< 5", "5-19", "20-49", "50-199", "≥ 200"]], index=True) | |
| 75 | +h2n = pd.DataFrame({k.split("|")[1]: v for k, v in H["H2"].items() if k.startswith("none")}).T; h2n.index.name = "densité"; h2n.columns = [f"k = {c}" for c in h2n.columns] | |
| 76 | +h2n_md = md_table(h2n.loc[["< 5", "5-19", "20-49", "50-199", "≥ 200"]], index=True) | |
| 77 | +h3 = pd.DataFrame({k.split("|")[3]: v for k, v in H["H3"].items() if k.startswith("5|5|12")}).T.loc[["< 5", "5-19", "20-49", "50-199", "≥ 200"]] | |
| 78 | +h3.columns = [LAB[c] for c in h3.columns]; h3.index.name = "densité"; h3_md = md_table(h3, index=True) | |
| 79 | +h3s = pd.DataFrame({k.split("|")[3]: v for k, v in H["H3"].items() if k.startswith("3|2|6")}).T.loc[["< 5", "5-19", "20-49", "50-199", "≥ 200"]] | |
| 80 | +h3s.columns = [LAB[c] for c in h3s.columns]; h3s.index.name = "densité"; h3s_md = md_table(h3s, index=True) | |
| 81 | +h4 = pd.DataFrame(H["H4"]).T.loc[["< 5", "5-19", "20-49", "50-199", "≥ 200"]]; h4.columns = [f"{LAB[c.split('|')[0]]} / {LAB[c.split('|')[1]]}" for c in h4.columns]; h4.index.name = "densité"; h4_md = md_table(h4, index=True) | |
| 82 | +h5 = pd.DataFrame(H["H5"]).T; h5.index.name = "sous-marché"; h5 = h5.rename(columns={"sd_median": "σ médian (log)"}); h5_md = md_table(h5[["σ médian (log)", "M2", "M7", "M8a", "M13", "M0", "M10", "M11b", "M12"]], index=True, fmt={"σ médian (log)": lambda v: f"{v:.3f}".replace(".", ",")}) | |
| 83 | +iv = pd.DataFrame(R["intervalles"]).T; iv.index = [f"{LAB[i.split('|')[0]]}, k = {i.split('|')[1]}, ≤ {i.split('|')[2]} km, {i.split('|')[3]} mois" for i in iv.index]; iv.index.name = "configuration" | |
| 84 | +iv_md = md_table(iv.rename(columns={"z80_val": "z calibré (val.)", "couverture_nominale_80": "couverture à z = 1,28", "couverture_calibree": "couverture calibrée", "largeur_mediane_pct": "largeur médiane (%)", "n_test": "n"}), index=True, fmt={"n": n_, "z calibré (val.)": f2}) | |
| 85 | +rc = pd.DataFrame(R["regles_conditionnelles"]).T.loc[["< 5", "5-19", "20-49", "50-199", "≥ 200"]] | |
| 86 | +rc["règle"] = rc.apply(lambda r: f"{LAB[r.sel]}, k = {int(r.k)}, ≤ {r.d} km, {int(r['T'])} mois, {LAB[r.adj]}, pond. {LAB[r.w]}, {LAB[r.recon]}", axis=1) | |
| 87 | +rc.index.name = "densité"; rc_md = md_table(rc.rename(columns={"mdape": "MdAPE val.", "mdape_test": "MdAPE test", "couverture": "couverture val.", "n_test": "n test"}), ["règle", "MdAPE val.", "MdAPE test", "couverture val.", "n test"], index=True, fmt={"n test": n_}) | |
| 88 | +stab = pd.DataFrame(R["stabilite"]); stab_md = md_table(stab[["trimestre", "n", "M2", "M7", "M0", "M11b", "M12", "biais_M7", "biais_M11b"]].rename(columns={"biais_M7": "biais M7", "biais_M11b": "biais M11b"}), fmt={"n": n_}) | |
| 89 | + | |
| 90 | +# ------------------------------------------------------------------ chiffres clés pour le texte | |
| 91 | +best = min(comp, key=lambda m: comp[m]["mdape"]); best_comps = min([m for m in comp if short(m) in ("M1", "M2", "M3", "M4", "M5", "M6", "M7", "M8a", "M8b", "M13", "M0", "M9")], key=lambda m: comp[m]["mdape"]) | |
| 92 | +m0 = R["M0"]; m0u = R["M0_sans_contrainte"] | |
| 93 | +gain_adj = comp["M2 moyenne des 5 plus proches (temps seul)"]["mdape"] - comp["M7 grille d'ajustement hédonique"]["mdape"] | |
| 94 | +gain_sel = comp["M2 moyenne des 5 plus proches (temps seul)"]["mdape"] - comp["M3 sélection par similarité physique"]["mdape"] | |
| 95 | +k_curve = {int(r.k): r.mdape_test for r in kb.itertuples()} | |
| 96 | +k_best = min(k_curve, key=k_curve.get) | |
| 97 | +e1b = R["E1_best_val"] | |
| 98 | +pair = R["paires"] | |
| 99 | +def pr(a, b): return pair[f"{a} vs {b}"] | |
| 100 | +M7 = "M7 grille d'ajustement hédonique"; M2 = "M2 moyenne des 5 plus proches (temps seul)"; M6 = "M6 grille d'ajustement conventionnelle" | |
| 101 | +M8a = "M8a grille locale (régression sur le voisinage)"; M8b = "M8b ajustements non linéaires (GBM)"; M11b = "M11b gradient boosting (avec rôle)"; M11a = "M11a gradient boosting (sans rôle)" | |
| 102 | +M0 = "M0 meilleure configuration empirique (choisie en 2024, couverture ≥ 97 %)"; M12 = "M12 hybride comparables + GBM (poids appris en 2024)"; M13 = "M13 ratio au rôle (5 plus proches, même rôle)"; M9 = "M9 réconciliation optimisée (poids appris en 2024)"; M3 = "M3 sélection par similarité physique"; M1 = "M1 médiane des 3 plus proches (temps seul)"; M10 = "M10 hédonique globale directe" | |
| 103 | +dens = R["strates_dens"]; densn = R["strates_dens_n"] | |
| 104 | +h1 = H["H1"] | |
| 105 | +marg = R["E2_marginal"] | |
| 106 | +cond_g = R["regles_conditionnelles_global_test"] | |
| 107 | +nT = R["n_targets"] | |
| 108 | + | |
| 109 | +txt = f"""--- | |
| 110 | +title: "La méthode des comparables à l'épreuve des données : validation empirique à grande échelle sur {n_(D['propre'])} ventes résidentielles québécoises" | |
| 111 | +subtitle: "Sélection, nombre, ajustements, pondération et réconciliation des ventes comparables — rétro-évaluation point-dans-le-temps sur les transactions Vrai Prix 2021-2026" | |
| 112 | +author: "Simon-Pierre Boucher — version {VERSION}" | |
| 113 | +date: "{time.strftime('%d/%m/%Y')}" | |
| 114 | +lang: fr-CA | |
| 115 | +--- | |
| 116 | + | |
| 117 | +# Résumé | |
| 118 | + | |
| 119 | +**Question.** La méthode de comparaison des ventes, telle qu'enseignée et pratiquée en évaluation immobilière, produit-elle des estimations fiables, et lesquelles de ses règles (choix des comparables, nombre, rayon, ancienneté, ajustements, pondération, réconciliation) sont empiriquement justifiées ? | |
| 120 | + | |
| 121 | +**Données et protocole.** {n_(D['propre'])} ventes de résidences unifamiliales du Québec (registre foncier jumelé au rôle d'évaluation, Vrai Prix, {D['dates'][0]} → {D['dates'][1]}), après un nettoyage journalisé de {n_(D['brut'])} transactions brutes. Pour chacune des **{n_(nT['val'] + nT['test'])} ventes cibles** de 2024 (validation, {n_(nT['val'])}) et de 2025-2026 (test hors échantillon, {n_(nT['test'])}), le prix est masqué et la valeur est estimée avec les seules ventes conclues **avant** la date de la cible, un indice de prix ré-estimé chaque mois avec les seules ventes antérieures, et des grilles d'ajustement et modèles ré-estimés chaque trimestre. **2 256 configurations** de la méthode des comparables et treize modèles de référence sont appliqués aux mêmes cibles ; les différences sont testées par bootstrap par grappes (municipalité × trimestre). | |
| 122 | + | |
| 123 | +**Résultats principaux (test 2025-2026).** | |
| 124 | + | |
| 125 | +1. **Les ajustements sont le facteur dominant.** Prendre la moyenne des cinq ventes les plus proches, ajustées seulement pour le temps (M2), donne une erreur absolue médiane (MdAPE) de **{pct(comp[M2]['mdape'])}** ; une grille d'ajustement hédonique (M7) la ramène à **{pct(comp[M7]['mdape'])}** (−{f1(gain_adj)} points, IC 95 % de la différence de MdAPE [{f2(pr(M7, M2)['ci_mdape'][0])} ; {f2(pr(M7, M2)['ci_mdape'][1])}]). Une grille conventionnelle stylisée (M6) fait presque aussi bien ({pct(comp[M6]['mdape'])}) : ce qui compte, c'est d'ajuster, plus que la finesse de la grille. | |
| 126 | +2. **Choisir des comparables semblables vaut autant qu'ajuster.** Sans aucun ajustement physique, sélectionner les cinq ventes les plus **semblables** (aire, terrain, âge, lien physique, étages) plutôt que les cinq plus **proches** fait passer la MdAPE de {pct(comp[M2]['mdape'])} à **{pct(comp[M3]['mdape'])}** (H1 confirmée : un comparable un peu plus éloigné mais semblable vaut mieux qu'un voisin dissemblable). | |
| 127 | +3. **Le nombre optimal de comparables est bien supérieur à trois.** Avec la grille hédonique, la MdAPE décroît de {pct(k_curve[1])} (k = 1) à {pct(k_curve[3])} (k = 3), {pct(k_curve[5])} (k = 5), {pct(k_curve[10])} (k = 10) et atteint son minimum à **k = {k_best}** ({pct(k_curve[k_best])}) avant de remonter ({pct(k_curve[50])} à k = 50) ; l'optimum dépend de la densité du marché (H2 confirmée). | |
| 128 | +4. **Le rayon et la fenêtre temporelle doivent être larges.** La meilleure configuration de la surface Erreur(k, d, t) choisie en 2024 est « {LAB[e1b['sel']]}, k = {int(e1b['k'])}, ≤ {e1b['d']} km, {int(e1b['T'])} mois, {LAB[e1b['adj']]} » ({pct(e1b['mdape_val'])} en validation → **{pct(e1b['mdape_test'])}** en test). Les rayons de 0,5-1 km et les fenêtres de 3 mois, privilégiés par la pratique, sont dominés dès que les comparables sont ajustés. | |
| 129 | +5. **La pondération et le mode de réconciliation comptent peu** (effets marginaux moyens de {f1(min(marg['w'].values()))} à {f1(max(marg['w'].values()))} points selon la pondération, médiane {f1(marg['recon']['median'])} contre moyenne {f1(marg['recon']['mean'])}), bien moins que la sélection ({f1(min(marg['sel'].values()))}-{f1(max(marg['sel'].values()))}) et l'ajustement ({f1(min(marg['adj'].values()))}-{f1(max(marg['adj'].values()))}). | |
| 130 | +6. **La meilleure méthode des comparables, entièrement transparente et choisie en 2024, atteint {pct(comp[M0]['mdape'])} en 2025-2026** (M0 : {LAB[m0['sel']]}, k = {int(m0['k'])}, ≤ {m0['d']} km, {int(m0['T'])} mois, {LAB[m0['adj']]}, pondération {LAB[m0['w']]}, {LAB[m0['recon']]}), contre **{pct(comp[M11b]['mdape'])}** pour un gradient boosting ré-entraîné chaque trimestre avec la valeur au rôle (M11b) et {pct(comp[M11a]['mdape'])} sans le rôle (M11a). L'hybride comparables + GBM (M12) atteint **{pct(comp[M12]['mdape'])}**. Le coût de la transparence est donc de {f1(comp[M0]['mdape'] - comp[M11b]['mdape'])} point{'s' if abs(comp[M0]['mdape'] - comp[M11b]['mdape']) >= 2 else ''} de MdAPE face au meilleur modèle opaque, et il s'annule dans les sous-marchés homogènes (H5). | |
| 131 | +7. **Le ratio au rôle est un comparable redoutable** : rôle de la cible × ratio prix/rôle médian des cinq plus proches ventes du même rôle (M13) donne {pct(comp[M13]['mdape'])}, mieux que toute grille d'ajustement physique, mais seulement pour {pct(comp[M13]['couverture'])} des cibles. | |
| 132 | +8. **La fiabilité dépend d'abord de la liquidité du sous-marché** : la MdAPE de M7 passe de {pct(dens['≥ 200']['M7'])} là où plus de 200 ventes ont eu lieu à moins de 2 km dans l'année à {pct(dens['< 5']['M7'])} là où il y en a eu moins de 5 ; des règles conditionnelles à la densité (système de recommandation) obtiennent {pct(cond_g['mdape'])} globalement avec {pct(cond_g['couverture'])} de couverture. | |
| 133 | + | |
| 134 | +**Portée.** Nous validons la précision prédictive de règles professionnelles *simulées* sur des caractéristiques observables (aire, terrain, âge, étages, lien physique, rôle) ; les ajustements non observables (garage, piscine, chambres, qualité, rénovations) restent dans l'erreur résiduelle, et le jugement réel d'un évaluateur n'est pas reproduit. | |
| 135 | + | |
| 136 | +# 1. Données | |
| 137 | + | |
| 138 | +## 1.1 Source et nettoyage | |
| 139 | + | |
| 140 | +Instantané de la base de production Vrai Prix pris le 2026-10-10 : {n_(766596)} transactions résidentielles et 3,75 M d'unités d'évaluation. Nous retenons les ventes de type « unifamilial » jumelées à une unité du rôle, et appliquons des règles de nettoyage explicites et journalisées — rien n'est éliminé silencieusement. | |
| 141 | + | |
| 142 | +{md_table(excl)} | |
| 143 | +Jeu propre : **{n_(D['propre'])} ventes** ({D['dates'][0]} → {D['dates'][1]}) ; prix médian {n_(D['prix']['50%'])} $, moyenne {n_(D['prix']['mean'])} $. | |
| 144 | + | |
| 145 | +Deux règles méritent un commentaire. R8 écarte les ratios prix/rôle extrêmes (< 0,33 ou > 3) : faute d'indicateur de lien de dépendance entre les parties, c'est notre proxy de « transaction non représentative » — un évaluateur les écarterait aussi. R6b écarte les ventes conclues l'année même de la construction : leur prix médian ({n_(346108)} $) et leur ratio au rôle (1,15) sont ceux de terrains ou de constructions inachevées, non de maisons achevées (450-480 k$ et 1,3 aux âges 1-4 ans). | |
| 146 | + | |
| 147 | +## 1.2 Répartition | |
| 148 | + | |
| 149 | +| Période | Rôle | Ventes | | |
| 150 | +|---|---|---| | |
| 151 | +| 2021-01 → 2023-12 | développement | {n_(D['par_periode']['train'])} | | |
| 152 | +| 2024 | validation (choix des configurations) | {n_(D['par_periode']['val'])} | | |
| 153 | +| 2025-01 → 2026-09 | test hors échantillon | {n_(D['par_periode']['test'])} | | |
| 154 | + | |
| 155 | +Régions les plus représentées : {', '.join(f"{k} ({n_(v)})" for k, v in list(D['par_region'].items())[:6])}. | |
| 156 | + | |
| 157 | +{{ width=85% }} | |
| 158 | + | |
| 159 | +## 1.3 Caractéristiques disponibles et limites | |
| 160 | + | |
| 161 | +Observées pour chaque vente : date, prix, coordonnées, municipalité, unité de voisinage, aire habitable, superficie du terrain, année de construction, nombre d'étages, lien physique (détaché, jumelé, en rangée), genre de construction, valeur au rôle en vigueur à la vente et sa part terrain. Non observées : garage, piscine, nombre de chambres et de salles de bains, qualité, état, rénovations. Les grilles d'ajustement testées ne portent donc que sur les premières ; les secondes sont une source d'erreur commune à toutes les méthodes. | |
| 162 | + | |
| 163 | +# 2. Protocole | |
| 164 | + | |
| 165 | +## 2.1 Rétro-évaluation point-dans-le-temps | |
| 166 | + | |
| 167 | +Pour chaque cible : (1) le prix est masqué ; (2) le bassin de comparables est constitué des ventes conclues strictement avant la date de la cible (la cible elle-même, si elle a déjà été vendue, est exclue : on évalue par comparables, non par vente répétée), dans un rayon de 25 km et une fenêtre de 24 mois ; (3) l'indice temporel régional utilisé pour ramener chaque comparable au niveau du dernier mois connu est ré-estimé chaque mois avec les seules ventes antérieures (régression hédonique à effets fixes de municipalité et effets mensuels, par région administrative, rétrécie vers l'indice provincial quand les ventes sont rares) ; (4) les coefficients d'ajustement hédoniques, les grilles régionales et les modèles d'apprentissage sont ré-estimés chaque trimestre avec les seules ventes antérieures au trimestre ; (5) chaque configuration produit une valeur et, le cas échéant, une dispersion des prix ajustés. | |
| 168 | + | |
| 169 | +Le niveau de prix du mois de la cible n'est jamais connu : toutes les méthodes travaillent au niveau du mois précédent, ce qui reproduit la situation réelle d'un évaluateur. | |
| 170 | + | |
| 171 | +## 2.2 Configurations testées | |
| 172 | + | |
| 173 | +**E1 — surface Erreur(k, d, t)** : k ∈ {{1, 3, 5, 7, 10, 15, 20, 30, 50}}, d ∈ {{0,5 ; 1 ; 2 ; 5 ; 10 ; 25 km}}, t ∈ {{3, 6, 12, 24 mois}}, sélection ∈ {{plus proches, plus semblables, compromis}}, ajustement ∈ {{temps seul, grille hédonique}}, médiane. 1 296 configurations. | |
| 174 | + | |
| 175 | +**E2 — plan factoriel des règles** aux ancres (k, d, t) = (3, 2 km, 6 mois), (5, 5 km, 12 mois), (10, 10 km, 24 mois) : sélection (4) × ajustement (8 : aucun, temps seul, grille conventionnelle, hédonique globale, hédonique régionale, régression locale, GBM, ratio au rôle) × pondération (5) × réconciliation (2). 960 configurations. | |
| 176 | + | |
| 177 | +*Grille conventionnelle stylisée* (M6) : aire habitable ajustée à la moitié du prix au m² du comparable ; terrain à la moitié de la valeur foncière implicite au m² (plafond ±15 %) ; âge à −0,5 % par année d'écart (plafond ±20 %) ; ±5 % entre détaché et attaché. *Grille hédonique* (M7) : coefficients d'une régression du log du prix sur les caractéristiques, à effets fixes de municipalité (within), ré-estimée chaque trimestre. *Régression locale* (M8a) : mêmes coefficients ré-estimés sur les ventes à moins de 10 km des 24 derniers mois (≥ 40 ventes, ridge vers la grille régionale). *GBM* (M8b) : l'ajustement est la différence des prédictions d'un gradient boosting entre la cible et le comparable. *Ratio au rôle* (M13) : valeur au rôle de la cible × ratio prix/rôle du comparable, pour les comparables de la même municipalité et du même rôle triennal. | |
| 178 | + | |
| 179 | +## 2.3 Mesures | |
| 180 | + | |
| 181 | +Erreur absolue en pourcentage APE = |P̂ − P| / P ; MdAPE (médiane), MAPE (moyenne), biais médian signé, parts à ±5, ±10 et ±20 %, couverture (part des cibles pour lesquelles la configuration produit une valeur — une configuration trop restrictive ne « rate » pas, elle s'abstient, et cela doit se voir). Les différences entre méthodes sont des différences d'erreurs appariées sur les mêmes cibles, avec intervalles de confiance à 95 % par bootstrap par grappes municipalité × trimestre (400 réplications). | |
| 182 | + | |
| 183 | +# 3. Résultats | |
| 184 | + | |
| 185 | +## 3.1 Combien de comparables ? À quelle distance ? De quelle ancienneté ? | |
| 186 | + | |
| 187 | +{{ width=100% }} | |
| 188 | + | |
| 189 | +{{ width=90% }} | |
| 190 | + | |
| 191 | +**Sélection par distance, grille hédonique, fenêtre 12 mois — MdAPE test (%)** | |
| 192 | + | |
| 193 | +{surf_tab('geo', 'hedo', 12)} | |
| 194 | +**Sélection par distance, temps seulement, fenêtre 12 mois — MdAPE test (%)** | |
| 195 | + | |
| 196 | +{surf_tab('geo', 'none', 12)} | |
| 197 | +**Sélection par similarité physique, grille hédonique, fenêtre 12 mois — MdAPE test (%)** | |
| 198 | + | |
| 199 | +{surf_tab('sim', 'hedo', 12)} | |
| 200 | +Lecture. (i) Sans ajustement, l'erreur décroît avec k et avec le rayon tant que l'on reste sous quelques kilomètres, puis se stabilise ; (ii) avec la grille hédonique, elle décroît nettement plus et l'optimum se déplace vers des k et des rayons plus grands ; (iii) la fenêtre de 3 mois est systématiquement dominée par 12-24 mois dès que l'on ajuste pour le temps. Meilleures configurations de E1 par sélection × ajustement (choisies en 2024, évaluées en 2025-2026) : | |
| 201 | + | |
| 202 | +{md_table(pd.DataFrame(R['E1_best_by_sel_adj']).T.reset_index().rename(columns={'index': 'sélection | ajustement', 'mdape_val': 'MdAPE val.', 'mdape_test': 'MdAPE test', 'couverture_test': 'couverture test'}), ['sélection | ajustement', 'k', 'd', 'T', 'MdAPE val.', 'MdAPE test', 'couverture test'], fmt={'k': lambda v: str(int(v)), 'T': lambda v: str(int(v)), 'd': lambda v: f1(v)})} | |
| 203 | +### k optimal selon la densité (H2) | |
| 204 | + | |
| 205 | +{h2_md} | |
| 206 | +*(sélection par distance, ≤ 10 km, 12 mois, grille hédonique ; MdAPE test en %, par nombre de ventes disponibles à moins de 2 km au cours des 12 mois précédents)* | |
| 207 | + | |
| 208 | +Sans ajustement : | |
| 209 | + | |
| 210 | +{h2n_md} | |
| 211 | +{{ width=100% }} | |
| 212 | + | |
| 213 | +## 3.2 Quels ajustements sont justifiés ? | |
| 214 | + | |
| 215 | +{{ width=90% }} | |
| 216 | + | |
| 217 | +{t_adj_md} | |
| 218 | +Trois enseignements. **L'ajustement temporel seul** vaut {f1(t_adj.set_index('adj').loc['raw', 'mdape_test'] - t_adj.set_index('adj').loc['none', 'mdape_test'])} point ; **n'importe quelle grille physique** vaut {f1(t_adj.set_index('adj').loc['none', 'mdape_test'] - t_adj.set_index('adj').loc['conv', 'mdape_test'])} à {f1(t_adj.set_index('adj').loc['none', 'mdape_test'] - t_adj.set_index('adj').loc['gbm', 'mdape_test'])} points de plus ; **les raffinements** (grille régionale, régression locale, ajustements non linéaires) apportent peu face à la grille hédonique globale, et la régression locale est même {('moins bonne' if t_adj.set_index('adj').loc['local', 'mdape_test'] > t_adj.set_index('adj').loc['hedo', 'mdape_test'] else 'meilleure')} ({pct(t_adj.set_index('adj').loc['local', 'mdape_test'])} contre {pct(t_adj.set_index('adj').loc['hedo', 'mdape_test'])}). | |
| 219 | + | |
| 220 | +### Grilles locales et petits échantillons (H3) | |
| 221 | + | |
| 222 | +MdAPE test (%) par type d'ajustement et densité, ancre (5, 5 km, 12 mois) : | |
| 223 | + | |
| 224 | +{h3_md} | |
| 225 | +Ancre (3, 2 km, 6 mois) : | |
| 226 | + | |
| 227 | +{h3s_md} | |
| 228 | +## 3.3 Comment choisir les comparables ? (H1) | |
| 229 | + | |
| 230 | +Ancre (5, 5 km, 12 mois), grille hédonique, médiane : | |
| 231 | + | |
| 232 | +{t_sel_md} | |
| 233 | +Test direct de H1 (k = 5, 12 mois, temps seulement) : les cinq voisins à moins de 1 km donnent une MdAPE de {pct(h1['geo_d1']['mdape_test'])} (couverture {pct(h1['geo_d1']['couverture_test'])}) ; les cinq ventes les plus semblables à moins de 10 km, {pct(h1['sim_d10']['mdape_test'])} (couverture {pct(h1['sim_d10']['couverture_test'])}) ; les cinq plus proches à moins de 10 km, {pct(h1['geo_d10']['mdape_test'])}. Avec la grille hédonique : plus proches {pct(h1['geo_d10_hedo']['mdape_test'])}, plus semblables {pct(h1['sim_d10_hedo']['mdape_test'])}, compromis {pct(h1['mix_d10_hedo']['mdape_test'])}. | |
| 234 | + | |
| 235 | +## 3.4 Pondération et réconciliation (H4) | |
| 236 | + | |
| 237 | +Ancre (5, 5 km, 12 mois), sélection par distance, grille hédonique : | |
| 238 | + | |
| 239 | +{t_wr_md} | |
| 240 | +Selon la densité (MdAPE test, %) : | |
| 241 | + | |
| 242 | +{h4_md} | |
| 243 | +Effets marginaux moyens (MdAPE test moyenne sur toutes les autres modalités de E2) : sélection {', '.join(f"{LAB[k]} {f1(v)}" for k, v in marg['sel'].items())} ; ajustement {', '.join(f"{LAB[k]} {f1(v)}" for k, v in marg['adj'].items())} ; pondération {', '.join(f"{LAB[k]} {f1(v)}" for k, v in marg['w'].items())} ; réconciliation {', '.join(f"{LAB[k]} {f1(v)}" for k, v in marg['recon'].items())}. | |
| 244 | + | |
| 245 | +## 3.5 Compétition de modèles | |
| 246 | + | |
| 247 | +{{ width=95% }} | |
| 248 | + | |
| 249 | +{t_comp} | |
| 250 | +M0 : {LAB[m0['sel']]}, k = {int(m0['k'])}, ≤ {m0['d']} km, {int(m0['T'])} mois, {LAB[m0['adj']]}, pondération {LAB[m0['w']]}, {LAB[m0['recon']]} (meilleure des 2 256 configurations en 2024 parmi celles qui couvrent ≥ 97 % des cibles ; sans cette contrainte, la meilleure est « {LAB[m0u['sel']]}, k = {int(m0u['k'])}, ≤ {m0u['d']} km, {int(m0u['T'])} mois, {LAB[m0u['adj']]} » à {pct(m0u['mdape_test'])} mais pour {pct(m0u['couverture_test'])} des cibles seulement). Poids appris de M9 : {', '.join(f"{LAB[k]} {f2(v)}" for k, v in R['M9_weights'].items() if v > 0)}. Poids de M12 : {', '.join(f"{k} {f2(v)}" for k, v in R['M12_weights'].items() if v > 0)}. | |
| 251 | + | |
| 252 | +### Différences appariées (test ; bootstrap par grappes municipalité × trimestre) | |
| 253 | + | |
| 254 | +{t_pairs} | |
| 255 | +Un Δ négatif signifie que A est plus précis que B ; l'intervalle qui exclut zéro signale une différence robuste aux dépendances spatiales et temporelles. | |
| 256 | + | |
| 257 | +## 3.6 Où la méthode tient-elle, où défaille-t-elle ? | |
| 258 | + | |
| 259 | +{{ width=85% }} | |
| 260 | + | |
| 261 | +**Par densité** (ventes à moins de 2 km dans les 12 mois précédents) — MdAPE test (%) : | |
| 262 | + | |
| 263 | +{strata_md('dens', 'densité')} | |
| 264 | +**Par région administrative** : | |
| 265 | + | |
| 266 | +{strata_md('region_nom', 'région')} | |
| 267 | +**Par quintile de prix** : | |
| 268 | + | |
| 269 | +{strata_md('prix_q', 'quintile de prix')} | |
| 270 | +**Par âge du bâtiment** : | |
| 271 | + | |
| 272 | +{strata_md('age_cl', 'âge')} | |
| 273 | +**Par lien physique** (0 détaché, 1 jumelé, 2 en rangée, 3 intégré) : | |
| 274 | + | |
| 275 | +{strata_md('lien', 'lien physique')} | |
| 276 | +### Marchés homogènes et hétérogènes (H5) | |
| 277 | + | |
| 278 | +Terciles de la dispersion (écart-type du log des prix ajustés) des dix ventes les plus proches à ≤ 10 km / 24 mois — MdAPE test (%) : | |
| 279 | + | |
| 280 | +{h5_md} | |
| 281 | +### Stabilité dans le temps | |
| 282 | + | |
| 283 | +{stab_md} | |
| 284 | +Le biais trimestriel commun à toutes les méthodes reflète le fait que le niveau de prix du mois courant est inconnu (toutes travaillent au niveau du mois précédent) et la saisonnalité des prix. | |
| 285 | + | |
| 286 | +## 3.7 Intervalles de valeur | |
| 287 | + | |
| 288 | +La dispersion des prix ajustés des comparables retenus fournit un intervalle de valeur. Avec z = 1,28 (couverture nominale de 80 % si les erreurs étaient normales), puis avec un z calibré sur la validation 2024 : | |
| 289 | + | |
| 290 | +{iv_md} | |
| 291 | +## 3.8 Vers un système de recommandation des comparables | |
| 292 | + | |
| 293 | +Règle choisie en 2024 pour chaque strate de densité (parmi les 2 256 configurations couvrant ≥ 95 % des cibles de la strate), évaluée en 2025-2026 : | |
| 294 | + | |
| 295 | +{rc_md} | |
| 296 | +Globalement, ces règles conditionnelles atteignent une MdAPE de **{pct(cond_g['mdape'])}** ({pct(cond_g['p10'])} des estimations à ±10 %, couverture {pct(cond_g['couverture'])}), contre {pct(comp[M0]['mdape'])} pour la meilleure règle unique. | |
| 297 | + | |
| 298 | +# 4. Discussion | |
| 299 | + | |
| 300 | +**Ce que la méthode des comparables fait bien.** Correctement paramétrée — comparables semblables, nombre d'une dizaine, rayon et fenêtre larges, ajustements hédoniques ou même conventionnels — elle atteint une précision ({pct(comp[M0]['mdape'])}) que ne dépasse le meilleur modèle d'apprentissage ré-entraîné chaque trimestre que de {f1(comp[M0]['mdape'] - comp[M11b]['mdape'])} point{'s' if abs(comp[M0]['mdape'] - comp[M11b]['mdape']) >= 2 else ''}, tout en restant intégralement explicable : chaque valeur est la médiane de quelques ventes identifiées, ajustées par des coefficients lisibles. | |
| 301 | + | |
| 302 | +**Ce qui contredit la pratique.** Trois comparables, un rayon de quelques centaines de mètres et une fenêtre de trois mois — les conventions les plus répandues — sont dominés sur toute la surface Erreur(k, d, t) : {pct(surf[(surf.sel == 'geo') & (surf.adj == 'hedo') & (surf.k == 3) & (surf.d == 1) & (surf['T'] == 3)].mdape_test.iloc[0])} de MdAPE (k = 3, ≤ 1 km, 3 mois, grille hédonique, couverture {pct(surf[(surf.sel == 'geo') & (surf.adj == 'hedo') & (surf.k == 3) & (surf.d == 1) & (surf['T'] == 3)].couverture_test.iloc[0])}) contre {pct(e1b['mdape_test'])} pour la meilleure configuration. La rareté des comparables « parfaits » est un argument pour élargir la recherche et ajuster, non pour se restreindre. | |
| 303 | + | |
| 304 | +**Le rôle d'évaluation comme information.** Le ratio au rôle (M13) et le GBM avec rôle (M11b) dominent leurs équivalents sans rôle : la valeur au rôle encapsule des caractéristiques non observées ici (qualité, dépendances, état) et une appréciation locale que les modèles globaux ne captent pas. Les modèles globaux sans rôle (M10, M11a) présentent d'ailleurs un biais négatif en 2025-2026 : l'appréciation est spatialement hétérogène et un indice régional, même ré-estimé chaque trimestre, la sous-estime là où les ventes se concentrent. | |
| 305 | + | |
| 306 | +**Limites.** (1) Règles professionnelles simulées, non dossiers d'évaluateurs. (2) Caractéristiques limitées à celles du rôle et du registre. (3) Délai de publication des ventes ignoré. (4) Un seul type de bien (unifamilial) ; condos et plex restent à traiter. (5) Les configurations ont été choisies sur 2024 puis évaluées sur 2025-2026 : les écarts entre validation et test (tableaux ci-dessus) mesurent la stabilité des règles, et ils sont faibles. | |
| 307 | + | |
| 308 | +# 5. Reproductibilité | |
| 309 | + | |
| 310 | +Code : `pipeline/01_dataset.py` → `06_report.py` (Python 3.14, pandas, numpy, scipy, lightgbm, matplotlib) ; sorties versionnées `out/{VERSION}/` (configurations `configs_E1.csv`, `configs_E2.csv` ; estimations `est_E1.npy`, `est_E2.npy` ; cibles `targets.parquet` ; tables `tables/*.csv` ; `results.json`). Exécution : `./run_all.sh` sur le nœud de calcul (M2U64), `./make_pdf.sh` sur le poste de travail. | |
| 311 | +""" | |
| 312 | +(DOCS / "RAPPORT.md").write_text(txt, encoding="utf-8") | |
| 313 | +log(f"rapport écrit : {DOCS / 'RAPPORT.md'} ({len(txt.split()):,} mots)") | |
added
pipeline/common.py
+85 −0
@@ -0,0 +1,85 @@ | ||
| 1 | +# -*- coding: utf-8 -*- | |
| 2 | +"""Validation empirique de la méthode des comparables — fonctions communes. | |
| 3 | +Données : instantané Vrai Prix (transactions ⋈ units) du 2026-10-10. | |
| 4 | +Toutes les sorties sont versionnées sous out/<version>/. | |
| 5 | +""" | |
| 6 | +from __future__ import annotations | |
| 7 | +import json, math, os, sqlite3, time | |
| 8 | +from pathlib import Path | |
| 9 | +import numpy as np | |
| 10 | +import pandas as pd | |
| 11 | + | |
| 12 | +ROOT = Path(os.environ.get("CV_ROOT", Path(__file__).resolve().parents[1])) | |
| 13 | +DB = Path(os.environ.get("CV_DB", ROOT / "data" / "vraiprix-20261010.db")) | |
| 14 | +REF = ROOT / "ref" | |
| 15 | +VERSION = os.environ.get("CV_VERSION", time.strftime("%Y-%m-%d")) | |
| 16 | +OUT = ROOT / "out" / VERSION | |
| 17 | +OUT.mkdir(parents=True, exist_ok=True) | |
| 18 | +(OUT / "fig").mkdir(exist_ok=True) | |
| 19 | + | |
| 20 | +# Périodes (séparation chronologique) | |
| 21 | +TRAIN_END = "2023-12-31" # développement | |
| 22 | +VAL_START, VAL_END = "2024-01-01", "2024-12-31" # choix des configurations | |
| 23 | +TEST_START, TEST_END = "2025-01-01", "2026-09-29" # évaluation finale hors échantillon | |
| 24 | +MONTH0 = "2021-01" # premier mois des données | |
| 25 | + | |
| 26 | +def log(msg: str) -> None: | |
| 27 | + print(f"[{time.strftime('%H:%M:%S')}] {msg}", flush=True) | |
| 28 | + | |
| 29 | +def connect(path: Path = DB) -> sqlite3.Connection: | |
| 30 | + return sqlite3.connect(f"file:{path}?mode=ro", uri=True, timeout=60) | |
| 31 | + | |
| 32 | +def _json_default(o): | |
| 33 | + if isinstance(o, np.integer): return int(o) | |
| 34 | + if isinstance(o, np.floating): return None if np.isnan(o) else float(o) | |
| 35 | + if isinstance(o, np.bool_): return bool(o) | |
| 36 | + if isinstance(o, np.ndarray): return o.tolist() | |
| 37 | + if isinstance(o, pd.Timestamp): return o.strftime("%Y-%m-%d") | |
| 38 | + if isinstance(o, float) and math.isnan(o): return None | |
| 39 | + return str(o) | |
| 40 | + | |
| 41 | +def save_json(obj, name: str) -> Path: | |
| 42 | + p = OUT / name; p.parent.mkdir(parents=True, exist_ok=True) | |
| 43 | + with open(p, "w", encoding="utf-8") as f: | |
| 44 | + json.dump(obj, f, ensure_ascii=False, indent=2, default=_json_default) | |
| 45 | + return p | |
| 46 | + | |
| 47 | +def load_json(name: str): | |
| 48 | + with open(OUT / name, encoding="utf-8") as f: | |
| 49 | + return json.load(f) | |
| 50 | + | |
| 51 | +def save_parquet(df: pd.DataFrame, name: str) -> Path: | |
| 52 | + p = OUT / name; df.to_parquet(p, index=False); return p | |
| 53 | + | |
| 54 | +def load_parquet(name: str) -> pd.DataFrame: | |
| 55 | + return pd.read_parquet(OUT / name) | |
| 56 | + | |
| 57 | +def month_index(dates) -> np.ndarray: | |
| 58 | + """Mois entier depuis MONTH0 (2021-01 = 0).""" | |
| 59 | + d = pd.to_datetime(dates) | |
| 60 | + y0, m0 = map(int, MONTH0.split("-")) | |
| 61 | + return ((d.dt.year - y0) * 12 + (d.dt.month - m0)).to_numpy() | |
| 62 | + | |
| 63 | +def to_xy_km(lat, lng, lat0=46.8, lng0=-71.2): | |
| 64 | + """Projection équirectangulaire locale (km) : erreur < 1 % pour des distances < 50 km.""" | |
| 65 | + x = (np.asarray(lng, float) - lng0) * 111.32 * np.cos(np.radians(np.asarray(lat, float))) | |
| 66 | + y = (np.asarray(lat, float) - lat0) * 110.57 | |
| 67 | + return x, y | |
| 68 | + | |
| 69 | +def load_repertoire() -> pd.DataFrame: | |
| 70 | + mun = pd.read_csv(REF / "MUN.csv", dtype=str, encoding="utf-8-sig") | |
| 71 | + mun = mun.rename(columns=str.lower) | |
| 72 | + out = pd.DataFrame({ | |
| 73 | + "code_mun": mun["mcode"].str.zfill(5), | |
| 74 | + "mun_nom": mun["munnom"], | |
| 75 | + "region": mun["regadm"].str.extract(r"\((\d+)\)")[0], | |
| 76 | + "region_nom": mun["regadm"].str.replace(r"\s*\(\d+\)", "", regex=True), | |
| 77 | + "mrc": mun["mrc"].str.extract(r"\((\d+)\)")[0], | |
| 78 | + "mrc_nom": mun["mrc"].str.replace(r"\s*\(\d+\)", "", regex=True), | |
| 79 | + }) | |
| 80 | + return out.drop_duplicates("code_mun") | |
| 81 | + | |
| 82 | +# Codages catégoriels (ordre stable) | |
| 83 | +LIEN = {"Détaché": 0, "Jumelé": 1, "En rangée 1 côté": 2, "En rangée plus de 1 côté": 2, "Intégré": 3} | |
| 84 | +GENRE = {"De plain-pied": 0, "Unimodulaire": 0, "À étages entiers": 1, "À niveaux décalés": 2, "À étage mansardé": 3} | |
| 85 | +BTYPE = {"single-story": 0, "full-story": 1, "mansard": 2, "split-level": 3, "single-wide": 4} | |
added
run_all.sh
+14 −0
@@ -0,0 +1,14 @@ | ||
| 1 | +#!/bin/zsh | |
| 2 | +# Chaîne complète sur le nœud de calcul (M2U64). Usage : CV_VERSION=2026-10-10 ./run_all.sh [étape de départ] | |
| 3 | +set -e | |
| 4 | +cd "$(dirname "$0")" | |
| 5 | +PY=${PY:-$HOME/uqo-eval-terrains/env/bin/python} | |
| 6 | +export CV_VERSION=${CV_VERSION:-$(date +%F)} | |
| 7 | +START=${1:-1} | |
| 8 | +for s in 01_dataset 02_index 03_models 04_backtest 05_analyze 06_report; do | |
| 9 | + n=${s%%_*} | |
| 10 | + [ "$((10#$n))" -lt "$START" ] && continue | |
| 11 | + echo "=== $s ($(date +%T)) ===" | |
| 12 | + (cd pipeline && $PY $s.py) | |
| 13 | +done | |
| 14 | +echo "=== terminé $(date +%T) ===" | |
| 15 | ||