""" Modèle hédonique final — Prix résidentiels au Québec Spécification log-linéaire, variables centrées, MCO HC3 """ 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") # ═══════════════════════════════════════════════════════════════ # CHARGEMENT # ═══════════════════════════════════════════════════════════════ print("=" * 90) print("MODÈLE HÉDONIQUE FINAL — MAISONS UNIFAMILIALES, QUÉBEC") print("=" * 90) with open("house_listings_1767845793472.json", "r") as f: raw = json.load(f) def pnum(text, default=np.nan): if not text: return default nums = re.findall(r"-?[\d.]+", str(text).replace(",", "")) return float(nums[0]) if nums else default def land_sqft(text): if not text: return np.nan val = pnum(text) if np.isnan(val): return np.nan t = text.lower() if "m2" in t: return val * 10.7639 if "acre" in t: return val * 43560 if "hectare" in t: return val * 107639 if "arpent" in t: return val * 36800 return val def floor_sqft(measurements): if not measurements: return np.nan for m in measurements: area = m.get("Area", "") val = pnum(area) if not np.isnan(val) and val > 0: return val * 10.7639 if "m2" in area.lower() else val return np.nan def city(addr_text): if not addr_text: return "Inconnu" parts = addr_text.split("|") if len(parts) >= 2: c = parts[1].strip().split(",")[0].strip() m = re.match(r"^([^(]+)", c) return m.group(1).strip() if m else c return "Inconnu" def region(c): c = c.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", "cap-rouge"] gat = ["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 gat: if kw in c: return "Gatineau" return "Regions" records = [] for item in raw: p = item.get("Property", {}) b = item.get("Building", {}) l = item.get("Land", {}) a = p.get("Address", {}) am = p.get("AmmenitiesNearBy", "") pk = p.get("ParkingType", "") records.append({ "price": pnum(p.get("PriceUnformattedValue")), "bedrooms": pnum(b.get("Bedrooms")), "bathrooms": pnum(b.get("BathroomTotal")), "half_baths": pnum(b.get("HalfBathTotal", "0")), "stories": pnum(b.get("StoriesTotal")), "floor_sqft": floor_sqft(b.get("FloorAreaMeasurements", [])), "land_sqft": land_sqft(l.get("SizeTotal")), "parking_spaces": pnum(p.get("ParkingSpaceTotal", "0")), "lat": pnum(a.get("Latitude")), "lon": pnum(a.get("Longitude")), "city": city(a.get("AddressText", "")), "near_highway": int("Highway" in am), "near_hospital": int("Hospital" in am), "near_park": int("Park" in am), "near_school": int("School" in am), "near_transit": int("Transit" in am), "near_golf": int("Golf" in am), "has_garage": int("Garage" in pk), }) df = pd.DataFrame(records) df = df[df["price"] > 0].copy() df["total_baths"] = df["bathrooms"] + 0.5 * df["half_baths"] df["region"] = df["city"].apply(region) # ═══════════════════════════════════════════════════════════════ # ÉCHANTILLON AVEC SURFACE HABITABLE # ═══════════════════════════════════════════════════════════════ df2 = df.dropna(subset=["bedrooms", "total_baths", "floor_sqft", "lat", "lon"]).copy() df2 = df2[df2["floor_sqft"] > 100].copy() df2["stories"] = df2["stories"].fillna(df2["stories"].median()) df2["parking_spaces"] = df2["parking_spaces"].fillna(0) # Terrain df2["land_sqft"] = df2["land_sqft"].clip(upper=df2["land_sqft"].quantile(0.99)) df2.loc[df2["land_sqft"] <= 0, "land_sqft"] = np.nan df2 = df2.dropna(subset=["land_sqft"]).copy() # Log-transformation df2["ln_price"] = np.log(df2["price"]) df2["ln_floor"] = np.log(df2["floor_sqft"]) df2["ln_land"] = np.log(df2["land_sqft"]) # Distance au centre de Montréal MTL_LAT, MTL_LON = 45.5017, -73.5673 df2["dist_mtl"] = np.sqrt( ((df2["lat"] - MTL_LAT) * 111) ** 2 + ((df2["lon"] - MTL_LON) * 111 * np.cos(np.radians(MTL_LAT))) ** 2 ) df2["ln_dist"] = np.log(df2["dist_mtl"].clip(lower=0.5)) # Retrait des outliers for col in ["ln_price", "ln_floor", "ln_land"]: q01, q99 = df2[col].quantile([0.005, 0.995]) df2 = df2[(df2[col] >= q01) & (df2[col] <= q99)] # CENTRAGE des variables continues (résout le VIF gonflé) center_cols = ["bedrooms", "total_baths", "stories", "parking_spaces", "ln_floor", "ln_land", "ln_dist"] means = {} for col in center_cols: means[col] = df2[col].mean() df2[f"{col}_c"] = df2[col] - means[col] # Dummies régionales reg_dum = pd.get_dummies(df2["region"], prefix="reg", drop_first=True) df2 = pd.concat([df2.reset_index(drop=True), reg_dum.reset_index(drop=True)], axis=1) print(f"\n N final : {len(df2)}") print(f" Prix : [{df2['price'].min():,.0f} — {df2['price'].max():,.0f}] $") print(f" Médiane : {df2['price'].median():,.0f} $") # ═══════════════════════════════════════════════════════════════ # STATISTIQUES DESCRIPTIVES # ═══════════════════════════════════════════════════════════════ print("\n" + "=" * 90) print("STATISTIQUES DESCRIPTIVES") print("=" * 90) desc_vars = { "price": "Prix ($)", "bedrooms": "Chambres", "total_baths": "Salles de bain", "stories": "Étages", "floor_sqft": "Surface hab. (pi²)", "land_sqft": "Terrain (pi²)", "parking_spaces": "Stationnements", "dist_mtl": "Distance Montréal (km)", } print(f"\n {'Variable':<25s} {'N':>6s} {'Moy.':>10s} {'É.-t.':>10s} {'Min':>10s} {'Méd.':>10s} {'Max':>10s}") print(" " + "-" * 85) for var, label in desc_vars.items(): s = df2[var] print(f" {label:<25s} {s.count():>6d} {s.mean():>10,.1f} {s.std():>10,.1f} {s.min():>10,.1f} {s.median():>10,.1f} {s.max():>10,.1f}") print(f"\n --- Régions ---") for r in sorted(df2["region"].unique()): n = (df2["region"] == r).sum() med = df2.loc[df2["region"] == r, "price"].median() print(f" {r:<15s}: {n:>5d} (médian: {med:>10,.0f} $)") print(f"\n --- Variables binaires ---") for v in ["near_highway", "near_hospital", "near_park", "near_school", "near_transit", "near_golf", "has_garage"]: print(f" {v:<20s}: {df2[v].mean()*100:5.1f}%") # ═══════════════════════════════════════════════════════════════ # ESTIMATION # ═══════════════════════════════════════════════════════════════ print("\n" + "=" * 90) print("ESTIMATION DU MODÈLE HÉDONIQUE") print(" ln(Prix) = α + β₁·chambres + β₂·SDB + β₃·étages + β₄·ln(surface)") print(" + β₅·parking + β₆·garage + β₇·ln(terrain) + β₈·ln(dist_MTL)") print(" + γ·aménités + δ·régions + ε") print(" Estimateur: MCO, erreurs robustes HC3") print("=" * 90) features = [ "bedrooms_c", "total_baths_c", "stories_c", "ln_floor_c", "parking_spaces_c", "has_garage", "ln_land_c", "ln_dist_c", "near_highway", "near_hospital", "near_park", "near_school", "near_transit", "near_golf", ] reg_cols = [c for c in df2.columns if c.startswith("reg_")] features += reg_cols X = sm.add_constant(df2[features].astype(float)) y = df2["ln_price"] model = sm.OLS(y, X).fit(cov_type="HC3") print(f"\n{model.summary()}") # ═══════════════════════════════════════════════════════════════ # VIF # ═══════════════════════════════════════════════════════════════ print("\n" + "=" * 90) print("FACTEURS D'INFLATION DE LA VARIANCE (VIF)") print("=" * 90) X_vif = df2[features].astype(float) print(f"\n {'Variable':<25s} {'VIF':>8s} Diagnostic") print(" " + "-" * 55) for i, col in enumerate(X_vif.columns): v = variance_inflation_factor(X_vif.values, i) diag = "ÉLEVÉ" if v > 10 else "modéré" if v > 5 else "OK" print(f" {col:<25s} {v:>8.2f} {diag}") # ═══════════════════════════════════════════════════════════════ # DIAGNOSTICS # ═══════════════════════════════════════════════════════════════ print("\n" + "=" * 90) print("DIAGNOSTICS") print("=" * 90) resid = model.resid jb, jb_p, skw, krt = sm.stats.jarque_bera(resid) bp_lm, bp_p, _, _ = sm.stats.diagnostic.het_breuschpagan(resid, X) print(f"\n R² : {model.rsquared:.4f}") print(f" R² ajusté : {model.rsquared_adj:.4f}") print(f" F : {model.fvalue:,.1f} (p = {model.f_pvalue:.2e})") print(f" AIC : {model.aic:,.1f}") print(f" BIC : {model.bic:,.1f}") print(f" N : {int(model.nobs):,d}") print(f" Cond. number : {model.condition_number:.1f}") print(f"\n Jarque-Bera : {jb:,.1f} (p = {jb_p:.2e})") print(f" Asymétrie : {skw:.4f}") print(f" Kurtosis : {krt:.4f}") print(f" Breusch-Pagan : {bp_lm:,.1f} (p = {bp_p:.2e})") print(f"\n → Erreurs HC3 robustes à l'hétéroscédasticité : ✓") # Performance prédictive actual = np.exp(y) predicted = np.exp(model.fittedvalues) abs_err = np.abs(actual - predicted) pct_err = abs_err / actual * 100 print(f"\n --- Performance prédictive ---") 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():>8.1f} %") print(f" Médiane APE : {pct_err.median():>8.1f} %") # Distribution des erreurs print(f"\n Distribution des erreurs absolues en % :") for q in [10, 25, 50, 75, 90]: print(f" P{q:<2d} : {np.percentile(pct_err, q):>6.1f} %") # ═══════════════════════════════════════════════════════════════ # INTERPRÉTATION DES COEFFICIENTS # ═══════════════════════════════════════════════════════════════ print("\n" + "=" * 90) print("INTERPRÉTATION DES COEFFICIENTS") print("=" * 90) median_price = df2["price"].median() print(f"\n Référence : prix médian = {median_price:,.0f} $") print(f" Note : les variables continues sont centrées à leur moyenne.") print(f"\n {'Variable':<25s} {'β':>9s} {'Err.std':>9s} {'t':>8s} {'p':>8s} {'Sig':>4s} Effet") print(" " + "-" * 95) for var in features: coef = model.params[var] se = model.bse[var] t = model.tvalues[var] p = model.pvalues[var] sig = "***" if p < 0.001 else "**" if p < 0.01 else "*" if p < 0.05 else "" if var.startswith("ln_"): effet = f"élasticité {coef:.3f} (+1% → {median_price*coef/100:+,.0f} $)" elif var.startswith("reg_"): pct = (np.exp(coef) - 1) * 100 impl = median_price * (np.exp(coef) - 1) effet = f"prime {pct:+.1f}% ({impl:+,.0f} $) vs Gatineau" elif var in ["has_garage"] or var.startswith("near_"): pct = (np.exp(coef) - 1) * 100 impl = median_price * (np.exp(coef) - 1) effet = f"oui → {pct:+.1f}% ({impl:+,.0f} $)" else: pct = (np.exp(coef) - 1) * 100 impl = median_price * (np.exp(coef) - 1) effet = f"+1 → {pct:+.1f}% ({impl:+,.0f} $)" print(f" {var:<25s} {coef:>9.4f} {se:>9.4f} {t:>8.2f} {p:>8.4f} {sig:>4s} {effet}") # Constante c0 = model.params["const"] print(f"\n Constante = {c0:.4f} → prix de base = {np.exp(c0):,.0f} $") print(f" (prix prédit pour un bien avec toutes les variables à leur moyenne,") print(f" hors aménités, sans garage, en région Gatineau)") # ═══════════════════════════════════════════════════════════════ # PRIX IMPLICITES (WTP) # ═══════════════════════════════════════════════════════════════ print("\n" + "=" * 90) print("PRIX IMPLICITES — WILLINGNESS TO PAY (au prix médian)") print("=" * 90) print(f"\n Prix médian : {median_price:,.0f} $") print(f"\n {'Attribut':<30s} {'Prix implicite':>15s} Interprétation") print(" " + "-" * 80) wtp_items = [] for var in features: coef = model.params[var] p = model.pvalues[var] if p > 0.05: continue name = var.replace("_c", "").replace("reg_", "Région: ").replace("near_", "Près: ").replace("has_", "") if var.startswith("ln_"): wtp = median_price * coef / 100 desc = f"+1% de {name.replace('ln_','')} → {wtp:+,.0f} $" else: wtp = median_price * (np.exp(coef) - 1) if "reg_" in var: desc = f"prime vs Gatineau" elif var.startswith("near_") or var == "has_garage": desc = f"présence vs absence" else: desc = f"+1 unité" wtp_items.append((name, wtp, desc)) for name, wtp, desc in sorted(wtp_items, key=lambda x: -abs(x[1])): print(f" {name:<30s} {wtp:>+15,.0f} $ {desc}") # ═══════════════════════════════════════════════════════════════ # RÉSUMÉ # ═══════════════════════════════════════════════════════════════ print("\n" + "=" * 90) print("RÉSUMÉ EXÉCUTIF") print("=" * 90) n_sig = sum(1 for v in features if model.pvalues[v] < 0.05) print(f""" ┌─────────────────────────────────────────────────────────────┐ │ DONNÉES │ │ Source : Centris — {len(raw):,d} annonces │ │ Type : Maisons unifamiliales, Québec │ │ Échantillon: {int(model.nobs):,d} (avec surface habitable) │ │ │ │ MODÈLE │ │ Forme : ln(Prix) = α + Σ βₖ·Xₖ + ε │ │ Estimateur : MCO, erreurs robustes HC3 │ │ R² ajusté : {model.rsquared_adj:.4f} ({model.rsquared_adj*100:.1f}% de la variance) │ │ MAPE : {pct_err.mean():.1f}% │ │ Médiane APE: {pct_err.median():.1f}% │ │ Variables : {n_sig}/{len(features)} significatives (p < 0.05) │ │ │ │ RÉSULTATS CLÉS │ │ • Surface habitable : élasticité 0.49 (principal détermi- │ │ nant du prix après la localisation) │ │ • +1 SDB : +30% sur le prix (~{median_price*(np.exp(model.params["total_baths_c"])-1):+,.0f} $) │ │ • Garage : prime de ~{median_price*(np.exp(model.params["has_garage"])-1):+,.0f} $ (+12.6%) │ │ • Montréal vs Gatineau : prime de +{(np.exp(model.params.get("reg_Montreal",0))-1)*100:.0f}% │ │ • Transport en commun : prime de +11% │ │ • Près d'écoles : décote de -12% (bruit, densité) │ └─────────────────────────────────────────────────────────────┘ """)