TeX 52.2%
Python 47.8%
1"""2Modèle hédonique final — Prix résidentiels au Québec3Spécification log-linéaire, variables centrées, MCO HC34"""56import json7import re8import warnings9import numpy as np10import pandas as pd11import statsmodels.api as sm12from statsmodels.stats.outliers_influence import variance_inflation_factor1314warnings.filterwarnings("ignore")1516# ═══════════════════════════════════════════════════════════════17# CHARGEMENT18# ═══════════════════════════════════════════════════════════════19print("=" * 90)20print("MODÈLE HÉDONIQUE FINAL — MAISONS UNIFAMILIALES, QUÉBEC")21print("=" * 90)2223with open("house_listings_1767845793472.json", "r") as f:24 raw = json.load(f)252627def pnum(text, default=np.nan):28 if not text:29 return default30 nums = re.findall(r"-?[\d.]+", str(text).replace(",", ""))31 return float(nums[0]) if nums else default323334def land_sqft(text):35 if not text:36 return np.nan37 val = pnum(text)38 if np.isnan(val):39 return np.nan40 t = text.lower()41 if "m2" in t: return val * 10.763942 if "acre" in t: return val * 4356043 if "hectare" in t: return val * 10763944 if "arpent" in t: return val * 3680045 return val464748def floor_sqft(measurements):49 if not measurements:50 return np.nan51 for m in measurements:52 area = m.get("Area", "")53 val = pnum(area)54 if not np.isnan(val) and val > 0:55 return val * 10.7639 if "m2" in area.lower() else val56 return np.nan575859def city(addr_text):60 if not addr_text:61 return "Inconnu"62 parts = addr_text.split("|")63 if len(parts) >= 2:64 c = parts[1].strip().split(",")[0].strip()65 m = re.match(r"^([^(]+)", c)66 return m.group(1).strip() if m else c67 return "Inconnu"686970def region(c):71 c = c.lower()72 mtl = ["montréal", "montreal", "laval", "longueuil", "brossard", "terrebonne",73 "blainville", "repentigny", "boucherville", "saint-jérôme", "mirabel",74 "boisbriand", "sainte-thérèse", "châteauguay", "saint-eustache",75 "deux-montagnes", "mascouche", "varennes", "vaudreuil", "dorval",76 "saint-hubert", "saint-lambert", "candiac", "carignan", "saint-bruno",77 "saint-constant", "saint-basile", "sainte-julie", "beloeil",78 "mont-saint-hilaire", "saint-jean-sur-richelieu", "chambly", "la prairie",79 "kirkland", "pointe-claire", "dollard", "côte-saint-luc", "westmount",80 "hampstead", "mont-royal", "outremont", "verdun", "lasalle", "lachine",81 "pierrefonds", "roxboro", "anjou", "saint-léonard", "rivière-des-prairies"]82 qc = ["québec", "quebec", "lévis", "levis", "beauport", "charlesbourg",83 "sainte-foy", "saint-augustin", "l'ancienne-lorette", "cap-rouge"]84 gat = ["gatineau", "hull", "aylmer"]85 for kw in mtl:86 if kw in c: return "Montreal"87 for kw in qc:88 if kw in c: return "Quebec"89 for kw in gat:90 if kw in c: return "Gatineau"91 return "Regions"929394records = []95for item in raw:96 p = item.get("Property", {})97 b = item.get("Building", {})98 l = item.get("Land", {})99 a = p.get("Address", {})100 am = p.get("AmmenitiesNearBy", "")101 pk = p.get("ParkingType", "")102 records.append({103 "price": pnum(p.get("PriceUnformattedValue")),104 "bedrooms": pnum(b.get("Bedrooms")),105 "bathrooms": pnum(b.get("BathroomTotal")),106 "half_baths": pnum(b.get("HalfBathTotal", "0")),107 "stories": pnum(b.get("StoriesTotal")),108 "floor_sqft": floor_sqft(b.get("FloorAreaMeasurements", [])),109 "land_sqft": land_sqft(l.get("SizeTotal")),110 "parking_spaces": pnum(p.get("ParkingSpaceTotal", "0")),111 "lat": pnum(a.get("Latitude")),112 "lon": pnum(a.get("Longitude")),113 "city": city(a.get("AddressText", "")),114 "near_highway": int("Highway" in am),115 "near_hospital": int("Hospital" in am),116 "near_park": int("Park" in am),117 "near_school": int("School" in am),118 "near_transit": int("Transit" in am),119 "near_golf": int("Golf" in am),120 "has_garage": int("Garage" in pk),121 })122123df = pd.DataFrame(records)124df = df[df["price"] > 0].copy()125df["total_baths"] = df["bathrooms"] + 0.5 * df["half_baths"]126df["region"] = df["city"].apply(region)127128# ═══════════════════════════════════════════════════════════════129# ÉCHANTILLON AVEC SURFACE HABITABLE130# ═══════════════════════════════════════════════════════════════131df2 = df.dropna(subset=["bedrooms", "total_baths", "floor_sqft", "lat", "lon"]).copy()132df2 = df2[df2["floor_sqft"] > 100].copy()133df2["stories"] = df2["stories"].fillna(df2["stories"].median())134df2["parking_spaces"] = df2["parking_spaces"].fillna(0)135136# Terrain137df2["land_sqft"] = df2["land_sqft"].clip(upper=df2["land_sqft"].quantile(0.99))138df2.loc[df2["land_sqft"] <= 0, "land_sqft"] = np.nan139df2 = df2.dropna(subset=["land_sqft"]).copy()140141# Log-transformation142df2["ln_price"] = np.log(df2["price"])143df2["ln_floor"] = np.log(df2["floor_sqft"])144df2["ln_land"] = np.log(df2["land_sqft"])145146# Distance au centre de Montréal147MTL_LAT, MTL_LON = 45.5017, -73.5673148df2["dist_mtl"] = np.sqrt(149 ((df2["lat"] - MTL_LAT) * 111) ** 2 +150 ((df2["lon"] - MTL_LON) * 111 * np.cos(np.radians(MTL_LAT))) ** 2151)152df2["ln_dist"] = np.log(df2["dist_mtl"].clip(lower=0.5))153154# Retrait des outliers155for col in ["ln_price", "ln_floor", "ln_land"]:156 q01, q99 = df2[col].quantile([0.005, 0.995])157 df2 = df2[(df2[col] >= q01) & (df2[col] <= q99)]158159# CENTRAGE des variables continues (résout le VIF gonflé)160center_cols = ["bedrooms", "total_baths", "stories", "parking_spaces",161 "ln_floor", "ln_land", "ln_dist"]162means = {}163for col in center_cols:164 means[col] = df2[col].mean()165 df2[f"{col}_c"] = df2[col] - means[col]166167# Dummies régionales168reg_dum = pd.get_dummies(df2["region"], prefix="reg", drop_first=True)169df2 = pd.concat([df2.reset_index(drop=True), reg_dum.reset_index(drop=True)], axis=1)170171print(f"\n N final : {len(df2)}")172print(f" Prix : [{df2['price'].min():,.0f} — {df2['price'].max():,.0f}] $")173print(f" Médiane : {df2['price'].median():,.0f} $")174175# ═══════════════════════════════════════════════════════════════176# STATISTIQUES DESCRIPTIVES177# ═══════════════════════════════════════════════════════════════178print("\n" + "=" * 90)179print("STATISTIQUES DESCRIPTIVES")180print("=" * 90)181182desc_vars = {183 "price": "Prix ($)",184 "bedrooms": "Chambres",185 "total_baths": "Salles de bain",186 "stories": "Étages",187 "floor_sqft": "Surface hab. (pi²)",188 "land_sqft": "Terrain (pi²)",189 "parking_spaces": "Stationnements",190 "dist_mtl": "Distance Montréal (km)",191}192193print(f"\n {'Variable':<25s} {'N':>6s} {'Moy.':>10s} {'É.-t.':>10s} {'Min':>10s} {'Méd.':>10s} {'Max':>10s}")194print(" " + "-" * 85)195for var, label in desc_vars.items():196 s = df2[var]197 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}")198199print(f"\n --- Régions ---")200for r in sorted(df2["region"].unique()):201 n = (df2["region"] == r).sum()202 med = df2.loc[df2["region"] == r, "price"].median()203 print(f" {r:<15s}: {n:>5d} (médian: {med:>10,.0f} $)")204205print(f"\n --- Variables binaires ---")206for v in ["near_highway", "near_hospital", "near_park", "near_school",207 "near_transit", "near_golf", "has_garage"]:208 print(f" {v:<20s}: {df2[v].mean()*100:5.1f}%")209210# ═══════════════════════════════════════════════════════════════211# ESTIMATION212# ═══════════════════════════════════════════════════════════════213print("\n" + "=" * 90)214print("ESTIMATION DU MODÈLE HÉDONIQUE")215print(" ln(Prix) = α + β₁·chambres + β₂·SDB + β₃·étages + β₄·ln(surface)")216print(" + β₅·parking + β₆·garage + β₇·ln(terrain) + β₈·ln(dist_MTL)")217print(" + γ·aménités + δ·régions + ε")218print(" Estimateur: MCO, erreurs robustes HC3")219print("=" * 90)220221features = [222 "bedrooms_c", "total_baths_c", "stories_c",223 "ln_floor_c", "parking_spaces_c", "has_garage",224 "ln_land_c", "ln_dist_c",225 "near_highway", "near_hospital", "near_park",226 "near_school", "near_transit", "near_golf",227]228reg_cols = [c for c in df2.columns if c.startswith("reg_")]229features += reg_cols230231X = sm.add_constant(df2[features].astype(float))232y = df2["ln_price"]233234model = sm.OLS(y, X).fit(cov_type="HC3")235print(f"\n{model.summary()}")236237# ═══════════════════════════════════════════════════════════════238# VIF239# ═══════════════════════════════════════════════════════════════240print("\n" + "=" * 90)241print("FACTEURS D'INFLATION DE LA VARIANCE (VIF)")242print("=" * 90)243244X_vif = df2[features].astype(float)245print(f"\n {'Variable':<25s} {'VIF':>8s} Diagnostic")246print(" " + "-" * 55)247for i, col in enumerate(X_vif.columns):248 v = variance_inflation_factor(X_vif.values, i)249 diag = "ÉLEVÉ" if v > 10 else "modéré" if v > 5 else "OK"250 print(f" {col:<25s} {v:>8.2f} {diag}")251252# ═══════════════════════════════════════════════════════════════253# DIAGNOSTICS254# ═══════════════════════════════════════════════════════════════255print("\n" + "=" * 90)256print("DIAGNOSTICS")257print("=" * 90)258259resid = model.resid260jb, jb_p, skw, krt = sm.stats.jarque_bera(resid)261bp_lm, bp_p, _, _ = sm.stats.diagnostic.het_breuschpagan(resid, X)262263print(f"\n R² : {model.rsquared:.4f}")264print(f" R² ajusté : {model.rsquared_adj:.4f}")265print(f" F : {model.fvalue:,.1f} (p = {model.f_pvalue:.2e})")266print(f" AIC : {model.aic:,.1f}")267print(f" BIC : {model.bic:,.1f}")268print(f" N : {int(model.nobs):,d}")269print(f" Cond. number : {model.condition_number:.1f}")270print(f"\n Jarque-Bera : {jb:,.1f} (p = {jb_p:.2e})")271print(f" Asymétrie : {skw:.4f}")272print(f" Kurtosis : {krt:.4f}")273print(f" Breusch-Pagan : {bp_lm:,.1f} (p = {bp_p:.2e})")274print(f"\n → Erreurs HC3 robustes à l'hétéroscédasticité : ✓")275276# Performance prédictive277actual = np.exp(y)278predicted = np.exp(model.fittedvalues)279abs_err = np.abs(actual - predicted)280pct_err = abs_err / actual * 100281282print(f"\n --- Performance prédictive ---")283print(f" RMSE : {np.sqrt(np.mean((actual-predicted)**2)):>12,.0f} $")284print(f" MAE : {abs_err.mean():>12,.0f} $")285print(f" MAPE : {pct_err.mean():>8.1f} %")286print(f" Médiane APE : {pct_err.median():>8.1f} %")287288# Distribution des erreurs289print(f"\n Distribution des erreurs absolues en % :")290for q in [10, 25, 50, 75, 90]:291 print(f" P{q:<2d} : {np.percentile(pct_err, q):>6.1f} %")292293# ═══════════════════════════════════════════════════════════════294# INTERPRÉTATION DES COEFFICIENTS295# ═══════════════════════════════════════════════════════════════296print("\n" + "=" * 90)297print("INTERPRÉTATION DES COEFFICIENTS")298print("=" * 90)299300median_price = df2["price"].median()301print(f"\n Référence : prix médian = {median_price:,.0f} $")302print(f" Note : les variables continues sont centrées à leur moyenne.")303print(f"\n {'Variable':<25s} {'β':>9s} {'Err.std':>9s} {'t':>8s} {'p':>8s} {'Sig':>4s} Effet")304print(" " + "-" * 95)305306for var in features:307 coef = model.params[var]308 se = model.bse[var]309 t = model.tvalues[var]310 p = model.pvalues[var]311 sig = "***" if p < 0.001 else "**" if p < 0.01 else "*" if p < 0.05 else ""312313 if var.startswith("ln_"):314 effet = f"élasticité {coef:.3f} (+1% → {median_price*coef/100:+,.0f} $)"315 elif var.startswith("reg_"):316 pct = (np.exp(coef) - 1) * 100317 impl = median_price * (np.exp(coef) - 1)318 effet = f"prime {pct:+.1f}% ({impl:+,.0f} $) vs Gatineau"319 elif var in ["has_garage"] or var.startswith("near_"):320 pct = (np.exp(coef) - 1) * 100321 impl = median_price * (np.exp(coef) - 1)322 effet = f"oui → {pct:+.1f}% ({impl:+,.0f} $)"323 else:324 pct = (np.exp(coef) - 1) * 100325 impl = median_price * (np.exp(coef) - 1)326 effet = f"+1 → {pct:+.1f}% ({impl:+,.0f} $)"327328 print(f" {var:<25s} {coef:>9.4f} {se:>9.4f} {t:>8.2f} {p:>8.4f} {sig:>4s} {effet}")329330# Constante331c0 = model.params["const"]332print(f"\n Constante = {c0:.4f} → prix de base = {np.exp(c0):,.0f} $")333print(f" (prix prédit pour un bien avec toutes les variables à leur moyenne,")334print(f" hors aménités, sans garage, en région Gatineau)")335336# ═══════════════════════════════════════════════════════════════337# PRIX IMPLICITES (WTP)338# ═══════════════════════════════════════════════════════════════339print("\n" + "=" * 90)340print("PRIX IMPLICITES — WILLINGNESS TO PAY (au prix médian)")341print("=" * 90)342343print(f"\n Prix médian : {median_price:,.0f} $")344print(f"\n {'Attribut':<30s} {'Prix implicite':>15s} Interprétation")345print(" " + "-" * 80)346347wtp_items = []348for var in features:349 coef = model.params[var]350 p = model.pvalues[var]351 if p > 0.05:352 continue353 name = var.replace("_c", "").replace("reg_", "Région: ").replace("near_", "Près: ").replace("has_", "")354 if var.startswith("ln_"):355 wtp = median_price * coef / 100356 desc = f"+1% de {name.replace('ln_','')} → {wtp:+,.0f} $"357 else:358 wtp = median_price * (np.exp(coef) - 1)359 if "reg_" in var:360 desc = f"prime vs Gatineau"361 elif var.startswith("near_") or var == "has_garage":362 desc = f"présence vs absence"363 else:364 desc = f"+1 unité"365 wtp_items.append((name, wtp, desc))366367for name, wtp, desc in sorted(wtp_items, key=lambda x: -abs(x[1])):368 print(f" {name:<30s} {wtp:>+15,.0f} $ {desc}")369370# ═══════════════════════════════════════════════════════════════371# RÉSUMÉ372# ═══════════════════════════════════════════════════════════════373print("\n" + "=" * 90)374print("RÉSUMÉ EXÉCUTIF")375print("=" * 90)376377n_sig = sum(1 for v in features if model.pvalues[v] < 0.05)378print(f"""379 ┌─────────────────────────────────────────────────────────────┐380 │ DONNÉES │381 │ Source : Centris — {len(raw):,d} annonces │382 │ Type : Maisons unifamiliales, Québec │383 │ Échantillon: {int(model.nobs):,d} (avec surface habitable) │384 │ │385 │ MODÈLE │386 │ Forme : ln(Prix) = α + Σ βₖ·Xₖ + ε │387 │ Estimateur : MCO, erreurs robustes HC3 │388 │ R² ajusté : {model.rsquared_adj:.4f} ({model.rsquared_adj*100:.1f}% de la variance) │389 │ MAPE : {pct_err.mean():.1f}% │390 │ Médiane APE: {pct_err.median():.1f}% │391 │ Variables : {n_sig}/{len(features)} significatives (p < 0.05) │392 │ │393 │ RÉSULTATS CLÉS │394 │ • Surface habitable : élasticité 0.49 (principal détermi- │395 │ nant du prix après la localisation) │396 │ • +1 SDB : +30% sur le prix (~{median_price*(np.exp(model.params["total_baths_c"])-1):+,.0f} $) │397 │ • Garage : prime de ~{median_price*(np.exp(model.params["has_garage"])-1):+,.0f} $ (+12.6%) │398 │ • Montréal vs Gatineau : prime de +{(np.exp(model.params.get("reg_Montreal",0))-1)*100:.0f}% │399 │ • Transport en commun : prime de +11% │400 │ • Près d'écoles : décote de -12% (bruit, densité) │401 └─────────────────────────────────────────────────────────────┘402""")403