""" Modèle hédonique des prix immobiliers résidentiels — Québec Source: 18 920 annonces Centris (maisons unifamiliales) Estimation par MCO (log-linéaire) avec diagnostics complets. """ import json import re import warnings import numpy as np import pandas as pd import statsmodels.api as sm from scipy import stats warnings.filterwarnings("ignore", category=FutureWarning) # ───────────────────────────────────────────────────────────────────── # 1. CHARGEMENT ET EXTRACTION DES VARIABLES # ───────────────────────────────────────────────────────────────────── print("=" * 80) print("MODÈLE HÉDONIQUE — PRIX DES MAISONS UNIFAMILIALES AU QUÉBEC") print("=" * 80) print("\n[1/6] Chargement des données...") with open("house_listings_1767845793472.json", "r") as f: raw = json.load(f) print(f" → {len(raw)} annonces chargées") 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 "acres" in t or "acre" in t: return val * 43560 elif "hectare" in t: return val * 107639 elif "arpent" in t: return val * 36800 return val # sqft par défaut 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): 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() # Nettoyer les parenthèses pour garder la ville principale match = re.match(r"^([^(]+)", city) if match: return match.group(1).strip() return city return "Inconnu" def extract_region(address_text): """Extraire le nom entre parenthèses (arrondissement / secteur).""" if not address_text: return "" match = re.search(r"\(([^)]+)\)", address_text) return match.group(1) if match else "" 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", "") rec = { "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), "postal_code": item.get("PostalCode", ""), # Amenities binaires "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, # Parking "has_garage": 1 if "Garage" in parking_type else 0, } records.append(rec) df = pd.DataFrame(records) print(f" → DataFrame: {df.shape[0]} lignes × {df.shape[1]} colonnes") # ───────────────────────────────────────────────────────────────────── # 2. NETTOYAGE ET STATISTIQUES DESCRIPTIVES # ───────────────────────────────────────────────────────────────────── print("\n[2/6] Nettoyage et statistiques descriptives...") # Filtrer les observations sans prix valide df = df[df["price"] > 0].copy() # Variable dépendante : log du prix df["ln_price"] = np.log(df["price"]) # Total bathrooms (full + half) df["total_baths"] = df["bathrooms"] + 0.5 * df["half_baths"] # Régions métropolitaines (basé sur les villes les plus fréquentes) top_cities = df["city"].value_counts().head(20) # Créer variable régionale : Montréal, Québec, Autres def classify_region(city): c = city.lower() montreal_keywords = [ "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" ] quebec_keywords = [ "québec", "quebec", "lévis", "levis", "beauport", "charlesbourg", "sainte-foy", "saint-augustin", "l'ancienne-lorette", "shannon", "stoneham", "saint-nicolas", "saint-romuald" ] for kw in montreal_keywords: if kw in c: return "RMR_Montreal" for kw in quebec_keywords: if kw in c: return "RMR_Quebec" return "Autres_Regions" df["region"] = df["city"].apply(classify_region) print(f"\n Observations après nettoyage: {len(df)}") print(f"\n --- Distribution des prix ---") print(f" Moyenne : {df['price'].mean():>12,.0f} $") print(f" Médiane : {df['price'].median():>12,.0f} $") print(f" Écart-type : {df['price'].std():>12,.0f} $") print(f" Min : {df['price'].min():>12,.0f} $") print(f" Max : {df['price'].max():>12,.0f} $") print(f" Q1 (25%) : {df['price'].quantile(0.25):>12,.0f} $") print(f" Q3 (75%) : {df['price'].quantile(0.75):>12,.0f} $") print(f"\n --- Distribution régionale ---") for reg, count in df["region"].value_counts().items(): med = df.loc[df["region"] == reg, "price"].median() print(f" {reg:<20s}: {count:>6d} listings (prix médian: {med:>10,.0f} $)") print(f"\n --- Variables continues ---") cont_vars = ["price", "bedrooms", "total_baths", "stories", "floor_sqft", "land_sqft", "parking_spaces"] desc = df[cont_vars].describe().T[["count", "mean", "std", "min", "50%", "max"]] desc.columns = ["N", "Moyenne", "Écart-type", "Min", "Médiane", "Max"] print(desc.to_string(float_format=lambda x: f"{x:,.1f}")) print(f"\n --- Variables binaires (proportion = 1) ---") bin_vars = ["near_highway", "near_hospital", "near_park", "near_school", "near_transit", "near_university", "near_golf", "has_garage"] for v in bin_vars: pct = df[v].mean() * 100 print(f" {v:<20s}: {pct:5.1f}%") # ───────────────────────────────────────────────────────────────────── # 3. MODÈLE HÉDONIQUE — SPÉCIFICATION LOG-LINÉAIRE # ───────────────────────────────────────────────────────────────────── print("\n" + "=" * 80) print("[3/6] ESTIMATION DU MODÈLE HÉDONIQUE (MCO log-linéaire)") print("=" * 80) # Préparer les variables model_df = df.dropna(subset=["ln_price", "bedrooms", "total_baths"]).copy() # Remplacer les NaN pour les variables optionnelles model_df["stories"] = model_df["stories"].fillna(model_df["stories"].median()) model_df["parking_spaces"] = model_df["parking_spaces"].fillna(0) model_df["has_floor_area"] = (~model_df["floor_sqft"].isna()).astype(int) model_df["floor_sqft_filled"] = model_df["floor_sqft"].fillna(0) model_df["has_land_size"] = (~model_df["land_sqft"].isna()).astype(int) model_df["land_sqft_filled"] = model_df["land_sqft"].fillna(0) # Log des surfaces (quand disponible) model_df["ln_floor_sqft"] = np.where( model_df["floor_sqft_filled"] > 0, np.log(model_df["floor_sqft_filled"]), 0 ) model_df["ln_land_sqft"] = np.where( model_df["land_sqft_filled"] > 0, np.log(model_df["land_sqft_filled"]), 0 ) # Dummies régionales region_dummies = pd.get_dummies(model_df["region"], prefix="region", drop_first=True) model_df = pd.concat([model_df, region_dummies], axis=1) # Supprimer les outliers extrêmes (< 1er percentile et > 99e percentile du prix) p01 = model_df["ln_price"].quantile(0.01) p99 = model_df["ln_price"].quantile(0.99) model_df = model_df[(model_df["ln_price"] >= p01) & (model_df["ln_price"] <= p99)].copy() print(f" Observations après retrait des outliers (1%/99%): {len(model_df)}") # Définir les régresseurs feature_cols = [ "bedrooms", "total_baths", "stories", "parking_spaces", "has_garage", "ln_floor_sqft", "has_floor_area", "ln_land_sqft", "has_land_size", "near_highway", "near_hospital", "near_park", "near_school", "near_transit", "near_university", "near_golf", ] # Ajouter les dummies régionales region_cols = [c for c in model_df.columns if c.startswith("region_")] feature_cols += region_cols X = model_df[feature_cols].astype(float) X = sm.add_constant(X) y = model_df["ln_price"] # Estimation MCO avec erreurs-types robustes (HC3) model = sm.OLS(y, X).fit(cov_type="HC3") print(f"\n{model.summary()}") # ───────────────────────────────────────────────────────────────────── # 4. INTERPRÉTATION DES COEFFICIENTS # ───────────────────────────────────────────────────────────────────── print("\n" + "=" * 80) print("[4/6] INTERPRÉTATION DES COEFFICIENTS (effet marginal en %)") print("=" * 80) print(f"\n {'Variable':<25s} {'Coeff':>8s} {'p-value':>10s} {'Impact %':>10s} Interprétation") print(" " + "-" * 80) for var in feature_cols: coef = model.params[var] pval = model.pvalues[var] sig = "" if pval < 0.001: sig = "***" elif pval < 0.01: sig = "**" elif pval < 0.05: sig = "*" # Pour les variables log, le coeff est une élasticité if var.startswith("ln_"): impact_str = f"élast. {coef:.3f}" interp = f"↑ 1% surface → ↑ {coef:.3f}% prix" else: # Pour les autres, (exp(β)-1)*100 donne l'effet en % pct = (np.exp(coef) - 1) * 100 impact_str = f"{pct:+.1f}%" if "region" in var: ref = [c for c in ["Autres_Regions", "RMR_Montreal", "RMR_Quebec"] if f"region_{c}" not in region_cols][0] interp = f"vs. {ref}" elif var in ["bedrooms", "total_baths", "stories", "parking_spaces"]: interp = f"+1 unité → {pct:+.1f}% sur le prix" else: interp = f"présence → {pct:+.1f}% sur le prix" print(f" {var:<25s} {coef:>8.4f} {pval:>8.4f}{sig:>3s} {impact_str:>10s} {interp}") print(f"\n Significativité: *** p<0.001, ** p<0.01, * p<0.05") # ───────────────────────────────────────────────────────────────────── # 5. DIAGNOSTICS DU MODÈLE # ───────────────────────────────────────────────────────────────────── print("\n" + "=" * 80) print("[5/6] DIAGNOSTICS DU MODÈLE") print("=" * 80) residuals = model.resid fitted = model.fittedvalues # R² et R² ajusté print(f"\n R² : {model.rsquared:.4f}") print(f" R² ajusté : {model.rsquared_adj:.4f}") print(f" F-statistique: {model.fvalue:.1f} (p = {model.f_pvalue:.2e})") print(f" N obs. : {int(model.nobs)}") print(f" K variables : {len(feature_cols)}") print(f" AIC : {model.aic:.1f}") print(f" BIC : {model.bic:.1f}") # Normalité des résidus jb_stat, jb_pval, skew, kurtosis = sm.stats.jarque_bera(residuals) print(f"\n --- Test de normalité des résidus (Jarque-Bera) ---") print(f" JB stat : {jb_stat:.1f}") print(f" p-value : {jb_pval:.4e}") print(f" Asymétrie : {skew:.4f}") print(f" Kurtosis : {kurtosis:.4f}") if jb_pval < 0.05: print(" → Résidus non-normaux (attendu avec N grand, erreurs HC3 corrigent)") # Breusch-Pagan test (hétéroscédasticité) bp_lm, bp_pval, bp_f, bp_fpval = sm.stats.diagnostic.het_breuschpagan(residuals, X) print(f"\n --- Test d'hétéroscédasticité (Breusch-Pagan) ---") print(f" LM stat : {bp_lm:.1f}") print(f" p-value : {bp_pval:.4e}") if bp_pval < 0.05: print(" → Hétéroscédasticité détectée (erreurs HC3 robustes utilisées)") else: print(" → Pas d'hétéroscédasticité significative") # Condition number (multicollinéarité) print(f"\n --- Multicollinéarité ---") print(f" Condition number: {model.condition_number:.0f}") if model.condition_number > 30: print(" → Attention: possible multicollinéarité (>30)") # RMSE et MAE en niveau de prix model_df_clean = model_df.loc[y.index].copy() model_df_clean["predicted_ln"] = fitted model_df_clean["predicted_price"] = np.exp(fitted) model_df_clean["actual_price"] = np.exp(y) model_df_clean["abs_error"] = np.abs(model_df_clean["actual_price"] - model_df_clean["predicted_price"]) model_df_clean["pct_error"] = model_df_clean["abs_error"] / model_df_clean["actual_price"] * 100 rmse = np.sqrt(np.mean((model_df_clean["actual_price"] - model_df_clean["predicted_price"]) ** 2)) mae = model_df_clean["abs_error"].mean() mape = model_df_clean["pct_error"].mean() median_ape = model_df_clean["pct_error"].median() print(f"\n --- Performance prédictive (en $ et %) ---") print(f" RMSE : {rmse:>12,.0f} $") print(f" MAE : {mae:>12,.0f} $") print(f" MAPE : {mape:>8.1f} %") print(f" Médiane APE : {median_ape:>8.1f} %") # ───────────────────────────────────────────────────────────────────── # 6. PRIX IMPLICITES (WILLINGNESS TO PAY) # ───────────────────────────────────────────────────────────────────── print("\n" + "=" * 80) print("[6/6] PRIX IMPLICITES (valeur marginale au prix médian)") print("=" * 80) median_price = df["price"].median() print(f"\n Prix médian de l'échantillon: {median_price:,.0f} $") print(f"\n {'Attribut':<25s} {'Prix implicite':>15s} Description") print(" " + "-" * 70) for var in feature_cols: coef = model.params[var] pval = model.pvalues[var] if pval > 0.05: continue if var.startswith("ln_"): # Élasticité: +1% de surface → coeff% de prix implicit = median_price * coef / 100 desc = f"+1% → {implicit:+,.0f} $ (élasticité)" elif var.startswith("region_"): pct = (np.exp(coef) - 1) * 100 implicit = median_price * (np.exp(coef) - 1) desc = f"prime régionale: {implicit:+,.0f} $" else: implicit = median_price * (np.exp(coef) - 1) desc = f"+1 unité → {implicit:+,.0f} $" print(f" {var:<25s} {implicit:>+15,.0f} $ {desc}") # ───────────────────────────────────────────────────────────────────── # RÉSUMÉ FINAL # ───────────────────────────────────────────────────────────────────── print("\n" + "=" * 80) print("RÉSUMÉ") print("=" * 80) print(f""" Modèle : ln(Prix) = β₀ + Σ βₖ Xₖ + ε Méthode: MCO avec erreurs-types robustes HC3 N observations : {int(model.nobs):,d} R² ajusté : {model.rsquared_adj:.4f} MAPE : {mape:.1f}% Médiane APE : {median_ape:.1f}% Variables significatives (p < 0.05): {sum(model.pvalues[1:] < 0.05)} / {len(feature_cols)} Le modèle explique {model.rsquared_adj*100:.1f}% de la variance des log-prix. L'erreur médiane de prédiction est de {median_ape:.1f}%. """)