TeX 52%
Python 47.9%
1"""2Modèle hédonique des prix immobiliers résidentiels — Québec3Source: 18 920 annonces Centris (maisons unifamiliales)45Estimation par MCO (log-linéaire) avec diagnostics complets.6"""78import json9import re10import warnings11import numpy as np12import pandas as pd13import statsmodels.api as sm14from scipy import stats1516warnings.filterwarnings("ignore", category=FutureWarning)1718# ─────────────────────────────────────────────────────────────────────19# 1. CHARGEMENT ET EXTRACTION DES VARIABLES20# ─────────────────────────────────────────────────────────────────────2122print("=" * 80)23print("MODÈLE HÉDONIQUE — PRIX DES MAISONS UNIFAMILIALES AU QUÉBEC")24print("=" * 80)2526print("\n[1/6] Chargement des données...")27with open("house_listings_1767845793472.json", "r") as f:28 raw = json.load(f)29print(f" → {len(raw)} annonces chargées")303132def parse_numeric(text, default=np.nan):33 if not text:34 return default35 nums = re.findall(r"-?[\d.]+", str(text).replace(",", ""))36 if nums:37 try:38 return float(nums[0])39 except ValueError:40 return default41 return default424344def parse_land_sqft(text):45 if not text:46 return np.nan47 val = parse_numeric(text)48 if np.isnan(val):49 return np.nan50 t = text.lower()51 if "m2" in t:52 return val * 10.763953 elif "acres" in t or "acre" in t:54 return val * 4356055 elif "hectare" in t:56 return val * 10763957 elif "arpent" in t:58 return val * 3680059 return val # sqft par défaut606162def parse_floor_sqft(measurements):63 if not measurements:64 return np.nan65 for m in measurements:66 area = m.get("Area", "")67 val = parse_numeric(area)68 if not np.isnan(val):69 if "m2" in area.lower():70 return val * 10.763971 return val72 return np.nan737475def extract_city(address_text):76 if not address_text:77 return "Inconnu"78 parts = address_text.split("|")79 if len(parts) >= 2:80 city_part = parts[1].strip()81 city = city_part.split(",")[0].strip()82 # Nettoyer les parenthèses pour garder la ville principale83 match = re.match(r"^([^(]+)", city)84 if match:85 return match.group(1).strip()86 return city87 return "Inconnu"888990def extract_region(address_text):91 """Extraire le nom entre parenthèses (arrondissement / secteur)."""92 if not address_text:93 return ""94 match = re.search(r"\(([^)]+)\)", address_text)95 return match.group(1) if match else ""969798records = []99for item in raw:100 prop = item.get("Property", {})101 bld = item.get("Building", {})102 land = item.get("Land", {})103 addr = prop.get("Address", {})104 address_text = addr.get("AddressText", "")105106 amenities = prop.get("AmmenitiesNearBy", "")107 parking_type = prop.get("ParkingType", "")108109 rec = {110 "price": parse_numeric(prop.get("PriceUnformattedValue")),111 "bedrooms": parse_numeric(bld.get("Bedrooms")),112 "bathrooms": parse_numeric(bld.get("BathroomTotal")),113 "half_baths": parse_numeric(bld.get("HalfBathTotal", "0")),114 "stories": parse_numeric(bld.get("StoriesTotal")),115 "floor_sqft": parse_floor_sqft(bld.get("FloorAreaMeasurements", [])),116 "land_sqft": parse_land_sqft(land.get("SizeTotal")),117 "parking_spaces": parse_numeric(prop.get("ParkingSpaceTotal", "0")),118 "latitude": parse_numeric(addr.get("Latitude")),119 "longitude": parse_numeric(addr.get("Longitude")),120 "city": extract_city(address_text),121 "postal_code": item.get("PostalCode", ""),122 # Amenities binaires123 "near_highway": 1 if "Highway" in amenities else 0,124 "near_hospital": 1 if "Hospital" in amenities else 0,125 "near_park": 1 if "Park" in amenities else 0,126 "near_school": 1 if "School" in amenities else 0,127 "near_transit": 1 if "Transit" in amenities else 0,128 "near_university": 1 if "University" in amenities or "CEGEP" in amenities else 0,129 "near_golf": 1 if "Golf" in amenities else 0,130 # Parking131 "has_garage": 1 if "Garage" in parking_type else 0,132 }133 records.append(rec)134135df = pd.DataFrame(records)136print(f" → DataFrame: {df.shape[0]} lignes × {df.shape[1]} colonnes")137138# ─────────────────────────────────────────────────────────────────────139# 2. NETTOYAGE ET STATISTIQUES DESCRIPTIVES140# ─────────────────────────────────────────────────────────────────────141print("\n[2/6] Nettoyage et statistiques descriptives...")142143# Filtrer les observations sans prix valide144df = df[df["price"] > 0].copy()145146# Variable dépendante : log du prix147df["ln_price"] = np.log(df["price"])148149# Total bathrooms (full + half)150df["total_baths"] = df["bathrooms"] + 0.5 * df["half_baths"]151152# Régions métropolitaines (basé sur les villes les plus fréquentes)153top_cities = df["city"].value_counts().head(20)154155# Créer variable régionale : Montréal, Québec, Autres156def classify_region(city):157 c = city.lower()158 montreal_keywords = [159 "montréal", "montreal", "laval", "longueuil", "brossard",160 "terrebonne", "blainville", "repentigny", "boucherville",161 "saint-jérôme", "mirabel", "boisbriand", "sainte-thérèse",162 "châteauguay", "saint-eustache", "deux-montagnes",163 "mascouche", "varennes", "vaudreuil", "dorval",164 "saint-hubert", "saint-lambert", "candiac", "carignan",165 "saint-bruno", "saint-constant", "saint-basile",166 "sainte-julie", "beloeil", "mont-saint-hilaire",167 "saint-jean-sur-richelieu", "chambly"168 ]169 quebec_keywords = [170 "québec", "quebec", "lévis", "levis", "beauport",171 "charlesbourg", "sainte-foy", "saint-augustin",172 "l'ancienne-lorette", "shannon", "stoneham",173 "saint-nicolas", "saint-romuald"174 ]175 for kw in montreal_keywords:176 if kw in c:177 return "RMR_Montreal"178 for kw in quebec_keywords:179 if kw in c:180 return "RMR_Quebec"181 return "Autres_Regions"182183df["region"] = df["city"].apply(classify_region)184185print(f"\n Observations après nettoyage: {len(df)}")186print(f"\n --- Distribution des prix ---")187print(f" Moyenne : {df['price'].mean():>12,.0f} $")188print(f" Médiane : {df['price'].median():>12,.0f} $")189print(f" Écart-type : {df['price'].std():>12,.0f} $")190print(f" Min : {df['price'].min():>12,.0f} $")191print(f" Max : {df['price'].max():>12,.0f} $")192print(f" Q1 (25%) : {df['price'].quantile(0.25):>12,.0f} $")193print(f" Q3 (75%) : {df['price'].quantile(0.75):>12,.0f} $")194195print(f"\n --- Distribution régionale ---")196for reg, count in df["region"].value_counts().items():197 med = df.loc[df["region"] == reg, "price"].median()198 print(f" {reg:<20s}: {count:>6d} listings (prix médian: {med:>10,.0f} $)")199200print(f"\n --- Variables continues ---")201cont_vars = ["price", "bedrooms", "total_baths", "stories", "floor_sqft",202 "land_sqft", "parking_spaces"]203desc = df[cont_vars].describe().T[["count", "mean", "std", "min", "50%", "max"]]204desc.columns = ["N", "Moyenne", "Écart-type", "Min", "Médiane", "Max"]205print(desc.to_string(float_format=lambda x: f"{x:,.1f}"))206207print(f"\n --- Variables binaires (proportion = 1) ---")208bin_vars = ["near_highway", "near_hospital", "near_park", "near_school",209 "near_transit", "near_university", "near_golf", "has_garage"]210for v in bin_vars:211 pct = df[v].mean() * 100212 print(f" {v:<20s}: {pct:5.1f}%")213214# ─────────────────────────────────────────────────────────────────────215# 3. MODÈLE HÉDONIQUE — SPÉCIFICATION LOG-LINÉAIRE216# ─────────────────────────────────────────────────────────────────────217print("\n" + "=" * 80)218print("[3/6] ESTIMATION DU MODÈLE HÉDONIQUE (MCO log-linéaire)")219print("=" * 80)220221# Préparer les variables222model_df = df.dropna(subset=["ln_price", "bedrooms", "total_baths"]).copy()223224# Remplacer les NaN pour les variables optionnelles225model_df["stories"] = model_df["stories"].fillna(model_df["stories"].median())226model_df["parking_spaces"] = model_df["parking_spaces"].fillna(0)227model_df["has_floor_area"] = (~model_df["floor_sqft"].isna()).astype(int)228model_df["floor_sqft_filled"] = model_df["floor_sqft"].fillna(0)229model_df["has_land_size"] = (~model_df["land_sqft"].isna()).astype(int)230model_df["land_sqft_filled"] = model_df["land_sqft"].fillna(0)231232# Log des surfaces (quand disponible)233model_df["ln_floor_sqft"] = np.where(234 model_df["floor_sqft_filled"] > 0,235 np.log(model_df["floor_sqft_filled"]),236 0237)238model_df["ln_land_sqft"] = np.where(239 model_df["land_sqft_filled"] > 0,240 np.log(model_df["land_sqft_filled"]),241 0242)243244# Dummies régionales245region_dummies = pd.get_dummies(model_df["region"], prefix="region", drop_first=True)246model_df = pd.concat([model_df, region_dummies], axis=1)247248# Supprimer les outliers extrêmes (< 1er percentile et > 99e percentile du prix)249p01 = model_df["ln_price"].quantile(0.01)250p99 = model_df["ln_price"].quantile(0.99)251model_df = model_df[(model_df["ln_price"] >= p01) & (model_df["ln_price"] <= p99)].copy()252print(f" Observations après retrait des outliers (1%/99%): {len(model_df)}")253254# Définir les régresseurs255feature_cols = [256 "bedrooms",257 "total_baths",258 "stories",259 "parking_spaces",260 "has_garage",261 "ln_floor_sqft",262 "has_floor_area",263 "ln_land_sqft",264 "has_land_size",265 "near_highway",266 "near_hospital",267 "near_park",268 "near_school",269 "near_transit",270 "near_university",271 "near_golf",272]273274# Ajouter les dummies régionales275region_cols = [c for c in model_df.columns if c.startswith("region_")]276feature_cols += region_cols277278X = model_df[feature_cols].astype(float)279X = sm.add_constant(X)280y = model_df["ln_price"]281282# Estimation MCO avec erreurs-types robustes (HC3)283model = sm.OLS(y, X).fit(cov_type="HC3")284285print(f"\n{model.summary()}")286287# ─────────────────────────────────────────────────────────────────────288# 4. INTERPRÉTATION DES COEFFICIENTS289# ─────────────────────────────────────────────────────────────────────290print("\n" + "=" * 80)291print("[4/6] INTERPRÉTATION DES COEFFICIENTS (effet marginal en %)")292print("=" * 80)293print(f"\n {'Variable':<25s} {'Coeff':>8s} {'p-value':>10s} {'Impact %':>10s} Interprétation")294print(" " + "-" * 80)295296for var in feature_cols:297 coef = model.params[var]298 pval = model.pvalues[var]299 sig = ""300 if pval < 0.001:301 sig = "***"302 elif pval < 0.01:303 sig = "**"304 elif pval < 0.05:305 sig = "*"306307 # Pour les variables log, le coeff est une élasticité308 if var.startswith("ln_"):309 impact_str = f"élast. {coef:.3f}"310 interp = f"↑ 1% surface → ↑ {coef:.3f}% prix"311 else:312 # Pour les autres, (exp(β)-1)*100 donne l'effet en %313 pct = (np.exp(coef) - 1) * 100314 impact_str = f"{pct:+.1f}%"315 if "region" in var:316 ref = [c for c in ["Autres_Regions", "RMR_Montreal", "RMR_Quebec"]317 if f"region_{c}" not in region_cols][0]318 interp = f"vs. {ref}"319 elif var in ["bedrooms", "total_baths", "stories", "parking_spaces"]:320 interp = f"+1 unité → {pct:+.1f}% sur le prix"321 else:322 interp = f"présence → {pct:+.1f}% sur le prix"323324 print(f" {var:<25s} {coef:>8.4f} {pval:>8.4f}{sig:>3s} {impact_str:>10s} {interp}")325326print(f"\n Significativité: *** p<0.001, ** p<0.01, * p<0.05")327328# ─────────────────────────────────────────────────────────────────────329# 5. DIAGNOSTICS DU MODÈLE330# ─────────────────────────────────────────────────────────────────────331print("\n" + "=" * 80)332print("[5/6] DIAGNOSTICS DU MODÈLE")333print("=" * 80)334335residuals = model.resid336fitted = model.fittedvalues337338# R² et R² ajusté339print(f"\n R² : {model.rsquared:.4f}")340print(f" R² ajusté : {model.rsquared_adj:.4f}")341print(f" F-statistique: {model.fvalue:.1f} (p = {model.f_pvalue:.2e})")342print(f" N obs. : {int(model.nobs)}")343print(f" K variables : {len(feature_cols)}")344print(f" AIC : {model.aic:.1f}")345print(f" BIC : {model.bic:.1f}")346347# Normalité des résidus348jb_stat, jb_pval, skew, kurtosis = sm.stats.jarque_bera(residuals)349print(f"\n --- Test de normalité des résidus (Jarque-Bera) ---")350print(f" JB stat : {jb_stat:.1f}")351print(f" p-value : {jb_pval:.4e}")352print(f" Asymétrie : {skew:.4f}")353print(f" Kurtosis : {kurtosis:.4f}")354if jb_pval < 0.05:355 print(" → Résidus non-normaux (attendu avec N grand, erreurs HC3 corrigent)")356357# Breusch-Pagan test (hétéroscédasticité)358bp_lm, bp_pval, bp_f, bp_fpval = sm.stats.diagnostic.het_breuschpagan(residuals, X)359print(f"\n --- Test d'hétéroscédasticité (Breusch-Pagan) ---")360print(f" LM stat : {bp_lm:.1f}")361print(f" p-value : {bp_pval:.4e}")362if bp_pval < 0.05:363 print(" → Hétéroscédasticité détectée (erreurs HC3 robustes utilisées)")364else:365 print(" → Pas d'hétéroscédasticité significative")366367# Condition number (multicollinéarité)368print(f"\n --- Multicollinéarité ---")369print(f" Condition number: {model.condition_number:.0f}")370if model.condition_number > 30:371 print(" → Attention: possible multicollinéarité (>30)")372373# RMSE et MAE en niveau de prix374model_df_clean = model_df.loc[y.index].copy()375model_df_clean["predicted_ln"] = fitted376model_df_clean["predicted_price"] = np.exp(fitted)377model_df_clean["actual_price"] = np.exp(y)378model_df_clean["abs_error"] = np.abs(model_df_clean["actual_price"] - model_df_clean["predicted_price"])379model_df_clean["pct_error"] = model_df_clean["abs_error"] / model_df_clean["actual_price"] * 100380381rmse = np.sqrt(np.mean((model_df_clean["actual_price"] - model_df_clean["predicted_price"]) ** 2))382mae = model_df_clean["abs_error"].mean()383mape = model_df_clean["pct_error"].mean()384median_ape = model_df_clean["pct_error"].median()385386print(f"\n --- Performance prédictive (en $ et %) ---")387print(f" RMSE : {rmse:>12,.0f} $")388print(f" MAE : {mae:>12,.0f} $")389print(f" MAPE : {mape:>8.1f} %")390print(f" Médiane APE : {median_ape:>8.1f} %")391392# ─────────────────────────────────────────────────────────────────────393# 6. PRIX IMPLICITES (WILLINGNESS TO PAY)394# ─────────────────────────────────────────────────────────────────────395print("\n" + "=" * 80)396print("[6/6] PRIX IMPLICITES (valeur marginale au prix médian)")397print("=" * 80)398399median_price = df["price"].median()400print(f"\n Prix médian de l'échantillon: {median_price:,.0f} $")401print(f"\n {'Attribut':<25s} {'Prix implicite':>15s} Description")402print(" " + "-" * 70)403404for var in feature_cols:405 coef = model.params[var]406 pval = model.pvalues[var]407 if pval > 0.05:408 continue409410 if var.startswith("ln_"):411 # Élasticité: +1% de surface → coeff% de prix412 implicit = median_price * coef / 100413 desc = f"+1% → {implicit:+,.0f} $ (élasticité)"414 elif var.startswith("region_"):415 pct = (np.exp(coef) - 1) * 100416 implicit = median_price * (np.exp(coef) - 1)417 desc = f"prime régionale: {implicit:+,.0f} $"418 else:419 implicit = median_price * (np.exp(coef) - 1)420 desc = f"+1 unité → {implicit:+,.0f} $"421422 print(f" {var:<25s} {implicit:>+15,.0f} $ {desc}")423424# ─────────────────────────────────────────────────────────────────────425# RÉSUMÉ FINAL426# ─────────────────────────────────────────────────────────────────────427print("\n" + "=" * 80)428print("RÉSUMÉ")429print("=" * 80)430print(f"""431 Modèle : ln(Prix) = β₀ + Σ βₖ Xₖ + ε432 Méthode: MCO avec erreurs-types robustes HC3433434 N observations : {int(model.nobs):,d}435 R² ajusté : {model.rsquared_adj:.4f}436 MAPE : {mape:.1f}%437 Médiane APE : {median_ape:.1f}%438439 Variables significatives (p < 0.05): {sum(model.pvalues[1:] < 0.05)} / {len(feature_cols)}440441 Le modèle explique {model.rsquared_adj*100:.1f}% de la variance des log-prix.442 L'erreur médiane de prédiction est de {median_ape:.1f}%.443""")444