TeX 52%
Python 47.9%
1"""2Modèle hédonique des prix immobiliers — Québec (v2)3Spécification améliorée avec correction de la sélection de surface,4variables quadratiques, et interactions spatiales.5"""67import json8import re9import warnings10import numpy as np11import pandas as pd12import statsmodels.api as sm13from statsmodels.stats.outliers_influence import variance_inflation_factor1415warnings.filterwarnings("ignore")1617# ─────────────────────────────────────────────────────────────────────18# 1. CHARGEMENT ET EXTRACTION19# ─────────────────────────────────────────────────────────────────────20print("=" * 90)21print("MODÈLE HÉDONIQUE v2 — PRIX DES MAISONS UNIFAMILIALES AU QUÉBEC")22print("=" * 90)2324with open("house_listings_1767845793472.json", "r") as f:25 raw = json.load(f)262728def parse_numeric(text, default=np.nan):29 if not text:30 return default31 nums = re.findall(r"-?[\d.]+", str(text).replace(",", ""))32 if nums:33 try:34 return float(nums[0])35 except ValueError:36 return default37 return default383940def parse_land_sqft(text):41 if not text:42 return np.nan43 val = parse_numeric(text)44 if np.isnan(val):45 return np.nan46 t = text.lower()47 if "m2" in t:48 return val * 10.763949 elif "acre" in t:50 return val * 4356051 elif "hectare" in t:52 return val * 10763953 elif "arpent" in t:54 return val * 3680055 return val565758def parse_floor_sqft(measurements):59 if not measurements:60 return np.nan61 for m in measurements:62 area = m.get("Area", "")63 val = parse_numeric(area)64 if not np.isnan(val) and val > 0:65 if "m2" in area.lower():66 return val * 10.763967 return val68 return np.nan697071def extract_city(address_text):72 if not address_text:73 return "Inconnu"74 parts = address_text.split("|")75 if len(parts) >= 2:76 city_part = parts[1].strip()77 city = city_part.split(",")[0].strip()78 match = re.match(r"^([^(]+)", city)79 if match:80 return match.group(1).strip()81 return city82 return "Inconnu"838485def classify_region(city):86 c = city.lower()87 mtl = ["montréal", "montreal", "laval", "longueuil", "brossard",88 "terrebonne", "blainville", "repentigny", "boucherville",89 "saint-jérôme", "mirabel", "boisbriand", "sainte-thérèse",90 "châteauguay", "saint-eustache", "deux-montagnes",91 "mascouche", "varennes", "vaudreuil", "dorval",92 "saint-hubert", "saint-lambert", "candiac", "carignan",93 "saint-bruno", "saint-constant", "saint-basile",94 "sainte-julie", "beloeil", "mont-saint-hilaire",95 "saint-jean-sur-richelieu", "chambly", "la prairie",96 "kirkland", "pointe-claire", "dollard", "côte-saint-luc",97 "westmount", "hampstead", "mont-royal", "outremont",98 "verdun", "lasalle", "lachine", "pierrefonds", "roxboro",99 "anjou", "saint-léonard", "rivière-des-prairies"]100 qc = ["québec", "quebec", "lévis", "levis", "beauport",101 "charlesbourg", "sainte-foy", "saint-augustin",102 "l'ancienne-lorette", "shannon", "stoneham",103 "saint-nicolas", "saint-romuald", "cap-rouge"]104 sherb = ["sherbrooke", "magog", "fleurimont"]105 trois = ["trois-rivières", "trois-rivieres", "shawinigan", "cap-de-la-madeleine"]106 gatineau = ["gatineau", "hull", "aylmer"]107 for kw in mtl:108 if kw in c:109 return "Montreal"110 for kw in qc:111 if kw in c:112 return "Quebec"113 for kw in sherb:114 if kw in c:115 return "Sherbrooke"116 for kw in trois:117 if kw in c:118 return "TroisRivieres"119 for kw in gatineau:120 if kw in c:121 return "Gatineau"122 return "Regions"123124125records = []126for item in raw:127 prop = item.get("Property", {})128 bld = item.get("Building", {})129 land = item.get("Land", {})130 addr = prop.get("Address", {})131 address_text = addr.get("AddressText", "")132 amenities = prop.get("AmmenitiesNearBy", "")133 parking_type = prop.get("ParkingType", "")134135 records.append({136 "price": parse_numeric(prop.get("PriceUnformattedValue")),137 "bedrooms": parse_numeric(bld.get("Bedrooms")),138 "bathrooms": parse_numeric(bld.get("BathroomTotal")),139 "half_baths": parse_numeric(bld.get("HalfBathTotal", "0")),140 "stories": parse_numeric(bld.get("StoriesTotal")),141 "floor_sqft": parse_floor_sqft(bld.get("FloorAreaMeasurements", [])),142 "land_sqft": parse_land_sqft(land.get("SizeTotal")),143 "parking_spaces": parse_numeric(prop.get("ParkingSpaceTotal", "0")),144 "latitude": parse_numeric(addr.get("Latitude")),145 "longitude": parse_numeric(addr.get("Longitude")),146 "city": extract_city(address_text),147 "near_highway": 1 if "Highway" in amenities else 0,148 "near_hospital": 1 if "Hospital" in amenities else 0,149 "near_park": 1 if "Park" in amenities else 0,150 "near_school": 1 if "School" in amenities else 0,151 "near_transit": 1 if "Transit" in amenities else 0,152 "near_university": 1 if "University" in amenities or "CEGEP" in amenities else 0,153 "near_golf": 1 if "Golf" in amenities else 0,154 "has_garage": 1 if "Garage" in parking_type else 0,155 })156157df = pd.DataFrame(records)158df = df[df["price"] > 0].copy()159print(f"\n {len(df)} annonces chargées")160161# ─────────────────────────────────────────────────────────────────────162# 2. FEATURE ENGINEERING163# ─────────────────────────────────────────────────────────────────────164print("\n[2] Feature engineering...")165166df["ln_price"] = np.log(df["price"])167df["total_baths"] = df["bathrooms"] + 0.5 * df["half_baths"]168df["region"] = df["city"].apply(classify_region)169170# ── MODÈLE A : Sans surface habitable (N ≈ 18 000) ──────────────────171# ── MODÈLE B : Avec surface habitable (N ≈ 6 500, plus précis) ──────172173# Préparer le sous-échantillon complet174df_full = df.dropna(subset=["bedrooms", "total_baths", "latitude", "longitude"]).copy()175df_full["stories"] = df_full["stories"].fillna(df_full["stories"].median())176df_full["parking_spaces"] = df_full["parking_spaces"].fillna(0)177178# Traitement du terrain : clipper les valeurs aberrantes et log179df_full["land_sqft"] = df_full["land_sqft"].clip(upper=df_full["land_sqft"].quantile(0.99))180df_full.loc[df_full["land_sqft"] <= 0, "land_sqft"] = np.nan181df_full["ln_land"] = np.log(df_full["land_sqft"])182183# Distance au centre-ville de Montréal (45.5017, -73.5673)184MTL_LAT, MTL_LON = 45.5017, -73.5673185df_full["dist_montreal_km"] = np.sqrt(186 ((df_full["latitude"] - MTL_LAT) * 111) ** 2 +187 ((df_full["longitude"] - MTL_LON) * 111 * np.cos(np.radians(MTL_LAT))) ** 2188)189df_full["ln_dist_mtl"] = np.log(df_full["dist_montreal_km"].clip(lower=0.5))190191# Variables quadratiques192df_full["bedrooms_sq"] = df_full["bedrooms"] ** 2193df_full["baths_sq"] = df_full["total_baths"] ** 2194195# Interaction : chambres × salles de bain196df_full["bed_x_bath"] = df_full["bedrooms"] * df_full["total_baths"]197198# Nombre d'aménités proches (score d'urbanité)199amenity_cols = ["near_highway", "near_hospital", "near_park", "near_school",200 "near_transit", "near_university", "near_golf"]201df_full["amenity_score"] = df_full[amenity_cols].sum(axis=1)202203# Outliers204p01 = df_full["ln_price"].quantile(0.005)205p99 = df_full["ln_price"].quantile(0.995)206df_full = df_full[(df_full["ln_price"] >= p01) & (df_full["ln_price"] <= p99)].copy()207208# Sous-échantillon avec surface habitable209df_floor = df_full[df_full["floor_sqft"].notna() & (df_full["floor_sqft"] > 100)].copy()210df_floor["ln_floor"] = np.log(df_floor["floor_sqft"])211df_floor["floor_per_bed"] = df_floor["floor_sqft"] / df_floor["bedrooms"].clip(lower=1)212df_floor["ln_floor_per_bed"] = np.log(df_floor["floor_per_bed"])213214print(f" Modèle A (sans surface) : N = {len(df_full)}")215print(f" Modèle B (avec surface) : N = {len(df_floor)}")216217# ─────────────────────────────────────────────────────────────────────218# 3. STATISTIQUES DESCRIPTIVES219# ─────────────────────────────────────────────────────────────────────220print("\n" + "=" * 90)221print("[3] STATISTIQUES DESCRIPTIVES")222print("=" * 90)223224print(f"\n --- Distribution des prix ---")225for label, sub in [("Tous", df_full), ("Avec surface", df_floor)]:226 print(f"\n [{label}] N = {len(sub)}")227 print(f" Moyenne : {sub['price'].mean():>12,.0f} $")228 print(f" Médiane : {sub['price'].median():>12,.0f} $")229 print(f" Éc.-type : {sub['price'].std():>12,.0f} $")230 print(f" IQR : [{sub['price'].quantile(.25):,.0f} — {sub['price'].quantile(.75):,.0f}]")231232print(f"\n --- Régions ---")233for reg in sorted(df_full["region"].unique()):234 n = (df_full["region"] == reg).sum()235 med = df_full.loc[df_full["region"] == reg, "price"].median()236 print(f" {reg:<15s}: {n:>5d} (médian: {med:>10,.0f} $)")237238# ─────────────────────────────────────────────────────────────────────239# 4. MODÈLE A — SANS SURFACE HABITABLE240# ─────────────────────────────────────────────────────────────────────241print("\n" + "=" * 90)242print("[4] MODÈLE A — SANS SURFACE HABITABLE")243print("=" * 90)244245region_dummies_A = pd.get_dummies(df_full["region"], prefix="reg", drop_first=True)246df_A = pd.concat([df_full.reset_index(drop=True), region_dummies_A.reset_index(drop=True)], axis=1)247248features_A = [249 "bedrooms", "bedrooms_sq",250 "total_baths", "baths_sq",251 "bed_x_bath",252 "stories",253 "parking_spaces", "has_garage",254 "ln_land",255 "ln_dist_mtl",256 "near_highway", "near_hospital", "near_park",257 "near_school", "near_transit", "near_golf",258]259reg_cols_A = [c for c in df_A.columns if c.startswith("reg_")]260features_A += reg_cols_A261262df_A_clean = df_A.dropna(subset=features_A + ["ln_price"]).copy()263X_A = sm.add_constant(df_A_clean[features_A].astype(float))264y_A = df_A_clean["ln_price"]265266model_A = sm.OLS(y_A, X_A).fit(cov_type="HC3")267print(f"\n{model_A.summary()}")268269# ─────────────────────────────────────────────────────────────────────270# 5. MODÈLE B — AVEC SURFACE HABITABLE271# ─────────────────────────────────────────────────────────────────────272print("\n" + "=" * 90)273print("[5] MODÈLE B — AVEC SURFACE HABITABLE (échantillon réduit)")274print("=" * 90)275276region_dummies_B = pd.get_dummies(df_floor["region"], prefix="reg", drop_first=True)277df_B = pd.concat([df_floor.reset_index(drop=True), region_dummies_B.reset_index(drop=True)], axis=1)278279features_B = [280 "bedrooms",281 "total_baths",282 "stories",283 "ln_floor",284 "parking_spaces", "has_garage",285 "ln_land",286 "ln_dist_mtl",287 "near_highway", "near_hospital", "near_park",288 "near_school", "near_transit", "near_golf",289]290reg_cols_B = [c for c in df_B.columns if c.startswith("reg_")]291features_B += reg_cols_B292293df_B_clean = df_B.dropna(subset=features_B + ["ln_price"]).copy()294X_B = sm.add_constant(df_B_clean[features_B].astype(float))295y_B = df_B_clean["ln_price"]296297model_B = sm.OLS(y_B, X_B).fit(cov_type="HC3")298print(f"\n{model_B.summary()}")299300# ─────────────────────────────────────────────────────────────────────301# 6. COMPARAISON ET INTERPRÉTATION302# ─────────────────────────────────────────────────────────────────────303print("\n" + "=" * 90)304print("[6] COMPARAISON DES MODÈLES")305print("=" * 90)306307for label, mod, X, y in [("A (sans surface)", model_A, X_A, y_A),308 ("B (avec surface)", model_B, X_B, y_B)]:309 pred_ln = mod.fittedvalues310 actual = np.exp(y)311 predicted = np.exp(pred_ln)312 abs_err = np.abs(actual - predicted)313 pct_err = abs_err / actual * 100314315 print(f"\n Modèle {label}")316 print(f" N : {int(mod.nobs):,d}")317 print(f" R² : {mod.rsquared:.4f}")318 print(f" R² ajusté : {mod.rsquared_adj:.4f}")319 print(f" AIC : {mod.aic:,.0f}")320 print(f" BIC : {mod.bic:,.0f}")321 print(f" RMSE ($) : {np.sqrt(np.mean((actual - predicted)**2)):>12,.0f}")322 print(f" MAE ($) : {abs_err.mean():>12,.0f}")323 print(f" MAPE : {pct_err.mean():.1f}%")324 print(f" Médiane APE : {pct_err.median():.1f}%")325326# ─────────────────────────────────────────────────────────────────────327# 7. INTERPRÉTATION DÉTAILLÉE DU MEILLEUR MODÈLE328# ─────────────────────────────────────────────────────────────────────329print("\n" + "=" * 90)330print("[7] INTERPRÉTATION DES PRIX IMPLICITES")331print("=" * 90)332333best = model_B334best_features = features_B335median_price = df_floor["price"].median()336median_vals = df_B_clean[features_B].median()337338print(f"\n Modèle retenu : B (avec surface habitable)")339print(f" Prix médian : {median_price:,.0f} $")340print(f"\n {'Variable':<25s} {'β':>8s} {'p':>8s} {'Sig':>3s} {'Interprétation'}")341print(" " + "-" * 85)342343for var in best_features:344 coef = best.params[var]345 pval = best.pvalues[var]346 sig = "***" if pval < 0.001 else "**" if pval < 0.01 else "*" if pval < 0.05 else ""347348 if var.startswith("ln_"):349 # Élasticité350 interp = f"Élasticité = {coef:.3f} (+1% → {coef:.3f}% prix, soit ~{median_price*coef/100:+,.0f} $)"351 elif var.startswith("reg_"):352 pct = (np.exp(coef) - 1) * 100353 impl = median_price * (np.exp(coef) - 1)354 interp = f"Prime vs Régions : {pct:+.1f}% ({impl:+,.0f} $)"355 else:356 pct = (np.exp(coef) - 1) * 100357 impl = median_price * (np.exp(coef) - 1)358 interp = f"+1 unité → {pct:+.1f}% ({impl:+,.0f} $)"359360 print(f" {var:<25s} {coef:>8.4f} {pval:>8.4f} {sig:>3s} {interp}")361362# ─────────────────────────────────────────────────────────────────────363# 8. VIF (Variance Inflation Factors) pour le modèle B364# ─────────────────────────────────────────────────────────────────────365print("\n" + "=" * 90)366print("[8] FACTEURS D'INFLATION DE LA VARIANCE (VIF) — Modèle B")367print("=" * 90)368369X_vif = df_B_clean[features_B].astype(float)370vif_data = []371for i, col in enumerate(X_vif.columns):372 vif_val = variance_inflation_factor(X_vif.values, i)373 vif_data.append((col, vif_val))374375print(f"\n {'Variable':<25s} {'VIF':>8s} Note")376print(" " + "-" * 55)377for col, vif_val in sorted(vif_data, key=lambda x: -x[1]):378 note = ""379 if vif_val > 10:380 note = "⚠ ÉLEVÉ"381 elif vif_val > 5:382 note = "modéré"383 print(f" {col:<25s} {vif_val:>8.2f} {note}")384385# ─────────────────────────────────────────────────────────────────────386# 9. DIAGNOSTICS RÉSIDUS MODÈLE B387# ─────────────────────────────────────────────────────────────────────388print("\n" + "=" * 90)389print("[9] DIAGNOSTICS DES RÉSIDUS — Modèle B")390print("=" * 90)391392resid = model_B.resid393fitted = model_B.fittedvalues394395jb_stat, jb_pval, skew, kurt = sm.stats.jarque_bera(resid)396bp_lm, bp_pval, _, _ = sm.stats.diagnostic.het_breuschpagan(resid, X_B)397398print(f"\n Jarque-Bera : {jb_stat:,.1f} (p = {jb_pval:.2e})")399print(f" Asymétrie : {skew:.4f}")400print(f" Kurtosis : {kurt:.4f}")401print(f" Breusch-Pagan LM : {bp_lm:,.1f} (p = {bp_pval:.2e})")402print(f" Condition number : {model_B.condition_number:.0f}")403404if bp_pval < 0.05:405 print(" → Hétéroscédasticité présente — erreurs HC3 robustes utilisées ✓")406if jb_pval < 0.05:407 print(" → Non-normalité des résidus — inférence asymptotique valide (N grand) ✓")408if model_B.condition_number < 30:409 print(" → Pas de multicollinéarité problématique ✓")410411# Distribution des résidus par quantile412print(f"\n Résidus standardisés (distribution) :")413std_resid = resid / resid.std()414for q in [1, 5, 25, 50, 75, 95, 99]:415 val = np.percentile(std_resid, q)416 print(f" P{q:<2d} : {val:>+6.3f}")417418# ─────────────────────────────────────────────────────────────────────419# RÉSUMÉ EXÉCUTIF420# ─────────────────────────────────────────────────────────────────────421print("\n" + "=" * 90)422print("RÉSUMÉ EXÉCUTIF")423print("=" * 90)424425pred_B = np.exp(model_B.fittedvalues)426actual_B = np.exp(y_B)427mape_B = (np.abs(actual_B - pred_B) / actual_B * 100).mean()428medape_B = (np.abs(actual_B - pred_B) / actual_B * 100).median()429430print(f"""431 DONNÉES432 -------433 Source : Centris (annonces actives, Québec)434 Échantillon : {len(raw):,d} maisons unifamiliales435 Couverture : Province de Québec436437 MODÈLE RETENU (B — avec surface habitable)438 -------------------------------------------439 Spécification : ln(Prix) = β₀ + Σ βₖ·Xₖ + ε440 Estimateur : MCO avec erreurs robustes HC3441 N : {int(model_B.nobs):,d}442 R² ajusté : {model_B.rsquared_adj:.4f} ({model_B.rsquared_adj*100:.1f}% de la variance expliquée)443 MAPE : {mape_B:.1f}%444 Médiane APE : {medape_B:.1f}%445446 PRINCIPAUX RÉSULTATS (prix implicites au prix médian de {median_price:,.0f} $)447 ─────────────────────────────────────────────────────────────────────448""")449450key_vars = [v for v in best_features if best.pvalues[v] < 0.05 and not v.startswith("reg_")]451for var in key_vars:452 coef = best.params[var]453 if var.startswith("ln_"):454 desc = f"+1% → {median_price * coef / 100:+,.0f} $"455 else:456 impl = median_price * (np.exp(coef) - 1)457 desc = f"+1 unité → {impl:+,.0f} $"458 print(f" {var:<22s}: {desc}")459460reg_vars = [v for v in best_features if v.startswith("reg_") and best.pvalues[v] < 0.05]461if reg_vars:462 print(f"\n PRIMES RÉGIONALES (vs. Régions périphériques)")463 for var in reg_vars:464 coef = best.params[var]465 impl = median_price * (np.exp(coef) - 1)466 print(f" {var.replace('reg_', ''):<22s}: {impl:+,.0f} $ ({(np.exp(coef)-1)*100:+.1f}%)")467468print()469