""" Modèle hédonique des prix immobiliers — Québec (v2) Spécification améliorée avec correction de la sélection de surface, variables quadratiques, et interactions spatiales. """ import json import re import warnings import numpy as np import pandas as pd import statsmodels.api as sm from statsmodels.stats.outliers_influence import variance_inflation_factor warnings.filterwarnings("ignore") # ───────────────────────────────────────────────────────────────────── # 1. CHARGEMENT ET EXTRACTION # ───────────────────────────────────────────────────────────────────── print("=" * 90) print("MODÈLE HÉDONIQUE v2 — PRIX DES MAISONS UNIFAMILIALES AU QUÉBEC") print("=" * 90) with open("house_listings_1767845793472.json", "r") as f: raw = json.load(f) def parse_numeric(text, default=np.nan): if not text: return default nums = re.findall(r"-?[\d.]+", str(text).replace(",", "")) if nums: try: return float(nums[0]) except ValueError: return default return default def parse_land_sqft(text): if not text: return np.nan val = parse_numeric(text) if np.isnan(val): return np.nan t = text.lower() if "m2" in t: return val * 10.7639 elif "acre" in t: return val * 43560 elif "hectare" in t: return val * 107639 elif "arpent" in t: return val * 36800 return val def parse_floor_sqft(measurements): if not measurements: return np.nan for m in measurements: area = m.get("Area", "") val = parse_numeric(area) if not np.isnan(val) and val > 0: if "m2" in area.lower(): return val * 10.7639 return val return np.nan def extract_city(address_text): if not address_text: return "Inconnu" parts = address_text.split("|") if len(parts) >= 2: city_part = parts[1].strip() city = city_part.split(",")[0].strip() match = re.match(r"^([^(]+)", city) if match: return match.group(1).strip() return city return "Inconnu" def classify_region(city): c = city.lower() mtl = ["montréal", "montreal", "laval", "longueuil", "brossard", "terrebonne", "blainville", "repentigny", "boucherville", "saint-jérôme", "mirabel", "boisbriand", "sainte-thérèse", "châteauguay", "saint-eustache", "deux-montagnes", "mascouche", "varennes", "vaudreuil", "dorval", "saint-hubert", "saint-lambert", "candiac", "carignan", "saint-bruno", "saint-constant", "saint-basile", "sainte-julie", "beloeil", "mont-saint-hilaire", "saint-jean-sur-richelieu", "chambly", "la prairie", "kirkland", "pointe-claire", "dollard", "côte-saint-luc", "westmount", "hampstead", "mont-royal", "outremont", "verdun", "lasalle", "lachine", "pierrefonds", "roxboro", "anjou", "saint-léonard", "rivière-des-prairies"] qc = ["québec", "quebec", "lévis", "levis", "beauport", "charlesbourg", "sainte-foy", "saint-augustin", "l'ancienne-lorette", "shannon", "stoneham", "saint-nicolas", "saint-romuald", "cap-rouge"] sherb = ["sherbrooke", "magog", "fleurimont"] trois = ["trois-rivières", "trois-rivieres", "shawinigan", "cap-de-la-madeleine"] gatineau = ["gatineau", "hull", "aylmer"] for kw in mtl: if kw in c: return "Montreal" for kw in qc: if kw in c: return "Quebec" for kw in sherb: if kw in c: return "Sherbrooke" for kw in trois: if kw in c: return "TroisRivieres" for kw in gatineau: if kw in c: return "Gatineau" return "Regions" records = [] for item in raw: prop = item.get("Property", {}) bld = item.get("Building", {}) land = item.get("Land", {}) addr = prop.get("Address", {}) address_text = addr.get("AddressText", "") amenities = prop.get("AmmenitiesNearBy", "") parking_type = prop.get("ParkingType", "") records.append({ "price": parse_numeric(prop.get("PriceUnformattedValue")), "bedrooms": parse_numeric(bld.get("Bedrooms")), "bathrooms": parse_numeric(bld.get("BathroomTotal")), "half_baths": parse_numeric(bld.get("HalfBathTotal", "0")), "stories": parse_numeric(bld.get("StoriesTotal")), "floor_sqft": parse_floor_sqft(bld.get("FloorAreaMeasurements", [])), "land_sqft": parse_land_sqft(land.get("SizeTotal")), "parking_spaces": parse_numeric(prop.get("ParkingSpaceTotal", "0")), "latitude": parse_numeric(addr.get("Latitude")), "longitude": parse_numeric(addr.get("Longitude")), "city": extract_city(address_text), "near_highway": 1 if "Highway" in amenities else 0, "near_hospital": 1 if "Hospital" in amenities else 0, "near_park": 1 if "Park" in amenities else 0, "near_school": 1 if "School" in amenities else 0, "near_transit": 1 if "Transit" in amenities else 0, "near_university": 1 if "University" in amenities or "CEGEP" in amenities else 0, "near_golf": 1 if "Golf" in amenities else 0, "has_garage": 1 if "Garage" in parking_type else 0, }) df = pd.DataFrame(records) df = df[df["price"] > 0].copy() print(f"\n {len(df)} annonces chargées") # ───────────────────────────────────────────────────────────────────── # 2. FEATURE ENGINEERING # ───────────────────────────────────────────────────────────────────── print("\n[2] Feature engineering...") df["ln_price"] = np.log(df["price"]) df["total_baths"] = df["bathrooms"] + 0.5 * df["half_baths"] df["region"] = df["city"].apply(classify_region) # ── MODÈLE A : Sans surface habitable (N ≈ 18 000) ────────────────── # ── MODÈLE B : Avec surface habitable (N ≈ 6 500, plus précis) ────── # Préparer le sous-échantillon complet df_full = df.dropna(subset=["bedrooms", "total_baths", "latitude", "longitude"]).copy() df_full["stories"] = df_full["stories"].fillna(df_full["stories"].median()) df_full["parking_spaces"] = df_full["parking_spaces"].fillna(0) # Traitement du terrain : clipper les valeurs aberrantes et log df_full["land_sqft"] = df_full["land_sqft"].clip(upper=df_full["land_sqft"].quantile(0.99)) df_full.loc[df_full["land_sqft"] <= 0, "land_sqft"] = np.nan df_full["ln_land"] = np.log(df_full["land_sqft"]) # Distance au centre-ville de Montréal (45.5017, -73.5673) MTL_LAT, MTL_LON = 45.5017, -73.5673 df_full["dist_montreal_km"] = np.sqrt( ((df_full["latitude"] - MTL_LAT) * 111) ** 2 + ((df_full["longitude"] - MTL_LON) * 111 * np.cos(np.radians(MTL_LAT))) ** 2 ) df_full["ln_dist_mtl"] = np.log(df_full["dist_montreal_km"].clip(lower=0.5)) # Variables quadratiques df_full["bedrooms_sq"] = df_full["bedrooms"] ** 2 df_full["baths_sq"] = df_full["total_baths"] ** 2 # Interaction : chambres × salles de bain df_full["bed_x_bath"] = df_full["bedrooms"] * df_full["total_baths"] # Nombre d'aménités proches (score d'urbanité) amenity_cols = ["near_highway", "near_hospital", "near_park", "near_school", "near_transit", "near_university", "near_golf"] df_full["amenity_score"] = df_full[amenity_cols].sum(axis=1) # Outliers p01 = df_full["ln_price"].quantile(0.005) p99 = df_full["ln_price"].quantile(0.995) df_full = df_full[(df_full["ln_price"] >= p01) & (df_full["ln_price"] <= p99)].copy() # Sous-échantillon avec surface habitable df_floor = df_full[df_full["floor_sqft"].notna() & (df_full["floor_sqft"] > 100)].copy() df_floor["ln_floor"] = np.log(df_floor["floor_sqft"]) df_floor["floor_per_bed"] = df_floor["floor_sqft"] / df_floor["bedrooms"].clip(lower=1) df_floor["ln_floor_per_bed"] = np.log(df_floor["floor_per_bed"]) print(f" Modèle A (sans surface) : N = {len(df_full)}") print(f" Modèle B (avec surface) : N = {len(df_floor)}") # ───────────────────────────────────────────────────────────────────── # 3. STATISTIQUES DESCRIPTIVES # ───────────────────────────────────────────────────────────────────── print("\n" + "=" * 90) print("[3] STATISTIQUES DESCRIPTIVES") print("=" * 90) print(f"\n --- Distribution des prix ---") for label, sub in [("Tous", df_full), ("Avec surface", df_floor)]: print(f"\n [{label}] N = {len(sub)}") print(f" Moyenne : {sub['price'].mean():>12,.0f} $") print(f" Médiane : {sub['price'].median():>12,.0f} $") print(f" Éc.-type : {sub['price'].std():>12,.0f} $") print(f" IQR : [{sub['price'].quantile(.25):,.0f} — {sub['price'].quantile(.75):,.0f}]") print(f"\n --- Régions ---") for reg in sorted(df_full["region"].unique()): n = (df_full["region"] == reg).sum() med = df_full.loc[df_full["region"] == reg, "price"].median() print(f" {reg:<15s}: {n:>5d} (médian: {med:>10,.0f} $)") # ───────────────────────────────────────────────────────────────────── # 4. MODÈLE A — SANS SURFACE HABITABLE # ───────────────────────────────────────────────────────────────────── print("\n" + "=" * 90) print("[4] MODÈLE A — SANS SURFACE HABITABLE") print("=" * 90) region_dummies_A = pd.get_dummies(df_full["region"], prefix="reg", drop_first=True) df_A = pd.concat([df_full.reset_index(drop=True), region_dummies_A.reset_index(drop=True)], axis=1) features_A = [ "bedrooms", "bedrooms_sq", "total_baths", "baths_sq", "bed_x_bath", "stories", "parking_spaces", "has_garage", "ln_land", "ln_dist_mtl", "near_highway", "near_hospital", "near_park", "near_school", "near_transit", "near_golf", ] reg_cols_A = [c for c in df_A.columns if c.startswith("reg_")] features_A += reg_cols_A df_A_clean = df_A.dropna(subset=features_A + ["ln_price"]).copy() X_A = sm.add_constant(df_A_clean[features_A].astype(float)) y_A = df_A_clean["ln_price"] model_A = sm.OLS(y_A, X_A).fit(cov_type="HC3") print(f"\n{model_A.summary()}") # ───────────────────────────────────────────────────────────────────── # 5. MODÈLE B — AVEC SURFACE HABITABLE # ───────────────────────────────────────────────────────────────────── print("\n" + "=" * 90) print("[5] MODÈLE B — AVEC SURFACE HABITABLE (échantillon réduit)") print("=" * 90) region_dummies_B = pd.get_dummies(df_floor["region"], prefix="reg", drop_first=True) df_B = pd.concat([df_floor.reset_index(drop=True), region_dummies_B.reset_index(drop=True)], axis=1) features_B = [ "bedrooms", "total_baths", "stories", "ln_floor", "parking_spaces", "has_garage", "ln_land", "ln_dist_mtl", "near_highway", "near_hospital", "near_park", "near_school", "near_transit", "near_golf", ] reg_cols_B = [c for c in df_B.columns if c.startswith("reg_")] features_B += reg_cols_B df_B_clean = df_B.dropna(subset=features_B + ["ln_price"]).copy() X_B = sm.add_constant(df_B_clean[features_B].astype(float)) y_B = df_B_clean["ln_price"] model_B = sm.OLS(y_B, X_B).fit(cov_type="HC3") print(f"\n{model_B.summary()}") # ───────────────────────────────────────────────────────────────────── # 6. COMPARAISON ET INTERPRÉTATION # ───────────────────────────────────────────────────────────────────── print("\n" + "=" * 90) print("[6] COMPARAISON DES MODÈLES") print("=" * 90) for label, mod, X, y in [("A (sans surface)", model_A, X_A, y_A), ("B (avec surface)", model_B, X_B, y_B)]: pred_ln = mod.fittedvalues actual = np.exp(y) predicted = np.exp(pred_ln) abs_err = np.abs(actual - predicted) pct_err = abs_err / actual * 100 print(f"\n Modèle {label}") print(f" N : {int(mod.nobs):,d}") print(f" R² : {mod.rsquared:.4f}") print(f" R² ajusté : {mod.rsquared_adj:.4f}") print(f" AIC : {mod.aic:,.0f}") print(f" BIC : {mod.bic:,.0f}") print(f" RMSE ($) : {np.sqrt(np.mean((actual - predicted)**2)):>12,.0f}") print(f" MAE ($) : {abs_err.mean():>12,.0f}") print(f" MAPE : {pct_err.mean():.1f}%") print(f" Médiane APE : {pct_err.median():.1f}%") # ───────────────────────────────────────────────────────────────────── # 7. INTERPRÉTATION DÉTAILLÉE DU MEILLEUR MODÈLE # ───────────────────────────────────────────────────────────────────── print("\n" + "=" * 90) print("[7] INTERPRÉTATION DES PRIX IMPLICITES") print("=" * 90) best = model_B best_features = features_B median_price = df_floor["price"].median() median_vals = df_B_clean[features_B].median() print(f"\n Modèle retenu : B (avec surface habitable)") print(f" Prix médian : {median_price:,.0f} $") print(f"\n {'Variable':<25s} {'β':>8s} {'p':>8s} {'Sig':>3s} {'Interprétation'}") print(" " + "-" * 85) for var in best_features: coef = best.params[var] pval = best.pvalues[var] sig = "***" if pval < 0.001 else "**" if pval < 0.01 else "*" if pval < 0.05 else "" if var.startswith("ln_"): # Élasticité interp = f"Élasticité = {coef:.3f} (+1% → {coef:.3f}% prix, soit ~{median_price*coef/100:+,.0f} $)" elif var.startswith("reg_"): pct = (np.exp(coef) - 1) * 100 impl = median_price * (np.exp(coef) - 1) interp = f"Prime vs Régions : {pct:+.1f}% ({impl:+,.0f} $)" else: pct = (np.exp(coef) - 1) * 100 impl = median_price * (np.exp(coef) - 1) interp = f"+1 unité → {pct:+.1f}% ({impl:+,.0f} $)" print(f" {var:<25s} {coef:>8.4f} {pval:>8.4f} {sig:>3s} {interp}") # ───────────────────────────────────────────────────────────────────── # 8. VIF (Variance Inflation Factors) pour le modèle B # ───────────────────────────────────────────────────────────────────── print("\n" + "=" * 90) print("[8] FACTEURS D'INFLATION DE LA VARIANCE (VIF) — Modèle B") print("=" * 90) X_vif = df_B_clean[features_B].astype(float) vif_data = [] for i, col in enumerate(X_vif.columns): vif_val = variance_inflation_factor(X_vif.values, i) vif_data.append((col, vif_val)) print(f"\n {'Variable':<25s} {'VIF':>8s} Note") print(" " + "-" * 55) for col, vif_val in sorted(vif_data, key=lambda x: -x[1]): note = "" if vif_val > 10: note = "⚠ ÉLEVÉ" elif vif_val > 5: note = "modéré" print(f" {col:<25s} {vif_val:>8.2f} {note}") # ───────────────────────────────────────────────────────────────────── # 9. DIAGNOSTICS RÉSIDUS MODÈLE B # ───────────────────────────────────────────────────────────────────── print("\n" + "=" * 90) print("[9] DIAGNOSTICS DES RÉSIDUS — Modèle B") print("=" * 90) resid = model_B.resid fitted = model_B.fittedvalues jb_stat, jb_pval, skew, kurt = sm.stats.jarque_bera(resid) bp_lm, bp_pval, _, _ = sm.stats.diagnostic.het_breuschpagan(resid, X_B) print(f"\n Jarque-Bera : {jb_stat:,.1f} (p = {jb_pval:.2e})") print(f" Asymétrie : {skew:.4f}") print(f" Kurtosis : {kurt:.4f}") print(f" Breusch-Pagan LM : {bp_lm:,.1f} (p = {bp_pval:.2e})") print(f" Condition number : {model_B.condition_number:.0f}") if bp_pval < 0.05: print(" → Hétéroscédasticité présente — erreurs HC3 robustes utilisées ✓") if jb_pval < 0.05: print(" → Non-normalité des résidus — inférence asymptotique valide (N grand) ✓") if model_B.condition_number < 30: print(" → Pas de multicollinéarité problématique ✓") # Distribution des résidus par quantile print(f"\n Résidus standardisés (distribution) :") std_resid = resid / resid.std() for q in [1, 5, 25, 50, 75, 95, 99]: val = np.percentile(std_resid, q) print(f" P{q:<2d} : {val:>+6.3f}") # ───────────────────────────────────────────────────────────────────── # RÉSUMÉ EXÉCUTIF # ───────────────────────────────────────────────────────────────────── print("\n" + "=" * 90) print("RÉSUMÉ EXÉCUTIF") print("=" * 90) pred_B = np.exp(model_B.fittedvalues) actual_B = np.exp(y_B) mape_B = (np.abs(actual_B - pred_B) / actual_B * 100).mean() medape_B = (np.abs(actual_B - pred_B) / actual_B * 100).median() print(f""" DONNÉES ------- Source : Centris (annonces actives, Québec) Échantillon : {len(raw):,d} maisons unifamiliales Couverture : Province de Québec MODÈLE RETENU (B — avec surface habitable) ------------------------------------------- Spécification : ln(Prix) = β₀ + Σ βₖ·Xₖ + ε Estimateur : MCO avec erreurs robustes HC3 N : {int(model_B.nobs):,d} R² ajusté : {model_B.rsquared_adj:.4f} ({model_B.rsquared_adj*100:.1f}% de la variance expliquée) MAPE : {mape_B:.1f}% Médiane APE : {medape_B:.1f}% PRINCIPAUX RÉSULTATS (prix implicites au prix médian de {median_price:,.0f} $) ───────────────────────────────────────────────────────────────────── """) key_vars = [v for v in best_features if best.pvalues[v] < 0.05 and not v.startswith("reg_")] for var in key_vars: coef = best.params[var] if var.startswith("ln_"): desc = f"+1% → {median_price * coef / 100:+,.0f} $" else: impl = median_price * (np.exp(coef) - 1) desc = f"+1 unité → {impl:+,.0f} $" print(f" {var:<22s}: {desc}") reg_vars = [v for v in best_features if v.startswith("reg_") and best.pvalues[v] < 0.05] if reg_vars: print(f"\n PRIMES RÉGIONALES (vs. Régions périphériques)") for var in reg_vars: coef = best.params[var] impl = median_price * (np.exp(coef) - 1) print(f" {var.replace('reg_', ''):<22s}: {impl:+,.0f} $ ({(np.exp(coef)-1)*100:+.1f}%)") print()