""" robustness.py - Robustness checks for hedonic pricing analysis. Implements a battery of specification tests following best practice in applied micro-economics: - Outlier treatment (winsorisation, trimming, log, median regression) - Sub-sample regressions by room type and price segment - Progressive location control inclusion - Host-type stratification - Summary formatting for publication tables All functions return structured dicts that can be passed directly to ``tables.robustness_table()`` and ``tables.regression_comparison_table()``. Author: Simon-Pierre Boucher Date: 2026-05-20 """ from __future__ import annotations import logging from typing import Any, Dict, List, Optional, Tuple import numpy as np import pandas as pd logger = logging.getLogger(__name__) # --------------------------------------------------------------------------- # Internal helpers # --------------------------------------------------------------------------- def _safe_ols( y: pd.Series, X: pd.DataFrame, add_constant: bool = True, ) -> Any: """Run OLS and return results; raises on failure.""" import statsmodels.api as sm if add_constant: X = sm.add_constant(X, has_constant="skip") model = sm.OLS(y, X, missing="drop") return model.fit(cov_type="HC1") def _result_to_dict(res: Any, label: str = "") -> Dict[str, Any]: """Convert a statsmodels result to a flat dict for comparison tables.""" out: Dict[str, Any] = { "N": int(res.nobs), "R2": round(res.rsquared, 4), "Adj_R2": round(res.rsquared_adj, 4), "AIC": round(res.aic, 1), "BIC": round(res.bic, 1), } for v in res.params.index: out[v] = { "coef": round(float(res.params[v]), 6), "se": round(float(res.bse[v]), 6), "pvalue": round(float(res.pvalues[v]), 6), } return out def _winsorize_series(s: pd.Series, limits: Tuple[float, float] = (0.01, 0.01)) -> pd.Series: """Winsorize a series at the given lower/upper percentile limits.""" from scipy.stats.mstats import winsorize as _winz return pd.Series(_winz(s.dropna(), limits=limits), index=s.dropna().index) def _trim_series(s: pd.Series, lower: float = 0.01, upper: float = 0.99) -> pd.Series: """Trim observations outside the specified quantile range.""" lo = s.quantile(lower) hi = s.quantile(upper) mask = (s >= lo) & (s <= hi) return s[mask] # --------------------------------------------------------------------------- # 1. Winsorization / outlier treatment # --------------------------------------------------------------------------- def robustness_winsorization( df: pd.DataFrame, y_col: str, x_cols: List[str], limits: Tuple[float, float] = (0.01, 0.01), trim_quantiles: Tuple[float, float] = (0.01, 0.99), ) -> Dict[str, Dict[str, Any]]: """ Compare regressions under different outlier treatments. Specifications: 1. Raw OLS 2. Winsorized dependent variable (1st / 99th) 3. Trimmed sample 4. Log-transformed dependent variable 5. Median (LAD) regression Parameters ---------- df : pd.DataFrame y_col : str Dependent variable (price). x_cols : list[str] Independent variables. limits : tuple Winsorization percentile limits. trim_quantiles : tuple Lower and upper quantile bounds for trimming. Returns ------- dict ``{spec_name: results_dict}`` """ import statsmodels.api as sm logger.info("Robustness: winsorization / outlier treatment checks.") results: Dict[str, Dict[str, Any]] = {} sub = df.dropna(subset=[y_col] + x_cols).copy() X = sub[x_cols] y = sub[y_col] # 1. Raw OLS try: res_raw = _safe_ols(y, X) results["Raw OLS"] = _result_to_dict(res_raw) except Exception as e: logger.error("Raw OLS failed: %s", e) # 2. Winsorized try: y_wins = _winsorize_series(y, limits=limits) # Align X with winsorized y index X_wins = X.loc[y_wins.index] res_wins = _safe_ols(y_wins, X_wins) results["Winsorized (1%/99%)"] = _result_to_dict(res_wins) except Exception as e: logger.error("Winsorized OLS failed: %s", e) # 3. Trimmed try: lo = y.quantile(trim_quantiles[0]) hi = y.quantile(trim_quantiles[1]) mask = (y >= lo) & (y <= hi) res_trim = _safe_ols(y[mask], X.loc[mask]) results["Trimmed (1%-99%)"] = _result_to_dict(res_trim) except Exception as e: logger.error("Trimmed OLS failed: %s", e) # 4. Log-transformed try: pos_mask = y > 0 y_log = np.log(y[pos_mask]) res_log = _safe_ols(y_log, X.loc[pos_mask]) results["Log(price) OLS"] = _result_to_dict(res_log) except Exception as e: logger.error("Log OLS failed: %s", e) # 5. Median regression (LAD) try: X_const = sm.add_constant(X, has_constant="skip") model_med = sm.QuantReg(y, X_const) res_med = model_med.fit(q=0.5) med_dict: Dict[str, Any] = { "N": int(res_med.nobs), "R2": np.nan, # pseudo R2 not directly comparable "Adj_R2": np.nan, } # Pseudo R2 for quantile regression try: med_dict["Pseudo_R2"] = round(float(1 - res_med.prsquared), 4) except Exception: pass for v in res_med.params.index: med_dict[v] = { "coef": round(float(res_med.params[v]), 6), "se": round(float(res_med.bse[v]), 6), "pvalue": round(float(res_med.pvalues[v]), 6), } results["Median Regression"] = med_dict except Exception as e: logger.error("Median regression failed: %s", e) logger.info("Winsorization robustness: %d/%d specs succeeded.", len(results), 5) return results # --------------------------------------------------------------------------- # 2. By room type # --------------------------------------------------------------------------- def robustness_by_room_type( df: pd.DataFrame, y_col: str, x_cols: List[str], room_col: str = "property_type", min_obs: int = 50, ) -> Dict[str, Dict[str, Any]]: """ Separate regressions for each room/property type. Parameters ---------- df : pd.DataFrame y_col : str x_cols : list[str] room_col : str Column identifying room/property type. min_obs : int Minimum observations to run a regression for a given type. Returns ------- dict ``{room_type: results_dict}`` """ logger.info("Robustness: sub-sample by room type (%s).", room_col) sub = df.dropna(subset=[y_col] + x_cols + [room_col]).copy() results: Dict[str, Dict[str, Any]] = {} # Full sample first try: res_full = _safe_ols(sub[y_col], sub[x_cols]) results["Full Sample"] = _result_to_dict(res_full) except Exception as e: logger.error("Full sample OLS failed: %s", e) for rtype, grp in sub.groupby(room_col, observed=True): if len(grp) < min_obs: logger.warning("Skipping %s (n=%d < %d).", rtype, len(grp), min_obs) continue # Remove x_cols that have zero variance in this sub-sample x_valid = [c for c in x_cols if grp[c].nunique() > 1] if not x_valid: logger.warning("Skipping %s: no valid regressors.", rtype) continue try: res = _safe_ols(grp[y_col], grp[x_valid]) results[str(rtype)] = _result_to_dict(res) except Exception as e: logger.error("OLS for %s failed: %s", rtype, e) logger.info("Room-type robustness: %d sub-samples estimated.", len(results)) return results # --------------------------------------------------------------------------- # 3. By price segment # --------------------------------------------------------------------------- def robustness_by_price_segment( df: pd.DataFrame, y_col: str, x_cols: List[str], segments: Optional[Dict[str, Tuple[float, float]]] = None, price_category_col: Optional[str] = "price_category", min_obs: int = 50, ) -> Dict[str, Dict[str, Any]]: """ Separate regressions for budget, moderate, premium, luxury segments. Parameters ---------- df : pd.DataFrame y_col : str x_cols : list[str] segments : dict, optional ``{segment_name: (lower_bound, upper_bound)}``. If None, quartile-based segments are used. price_category_col : str, optional If present in df, use this column instead of numeric cuts. min_obs : int Minimum observations per segment. Returns ------- dict """ logger.info("Robustness: sub-sample by price segment.") sub = df.dropna(subset=[y_col] + x_cols).copy() results: Dict[str, Dict[str, Any]] = {} # Full sample try: res_full = _safe_ols(sub[y_col], sub[x_cols]) results["Full Sample"] = _result_to_dict(res_full) except Exception as e: logger.error("Full sample OLS failed: %s", e) if price_category_col and price_category_col in sub.columns: for seg, grp in sub.groupby(price_category_col, observed=True): if len(grp) < min_obs: continue x_valid = [c for c in x_cols if grp[c].nunique() > 1] if not x_valid: continue try: res = _safe_ols(grp[y_col], grp[x_valid]) results[str(seg)] = _result_to_dict(res) except Exception as e: logger.error("Segment %s failed: %s", seg, e) else: if segments is None: q25, q50, q75 = sub[y_col].quantile([0.25, 0.50, 0.75]) segments = { "Budget (Q1)": (sub[y_col].min(), q25), "Moderate (Q2)": (q25, q50), "Premium (Q3)": (q50, q75), "Luxury (Q4)": (q75, sub[y_col].max()), } for seg_name, (lo, hi) in segments.items(): grp = sub[(sub[y_col] >= lo) & (sub[y_col] <= hi)] if len(grp) < min_obs: logger.warning("Skipping %s (n=%d).", seg_name, len(grp)) continue x_valid = [c for c in x_cols if grp[c].nunique() > 1] if not x_valid: continue try: res = _safe_ols(grp[y_col], grp[x_valid]) results[seg_name] = _result_to_dict(res) except Exception as e: logger.error("Segment %s failed: %s", seg_name, e) logger.info("Price-segment robustness: %d segments estimated.", len(results)) return results # --------------------------------------------------------------------------- # 4. Location controls # --------------------------------------------------------------------------- def robustness_location_controls( df: pd.DataFrame, y_col: str, x_cols: List[str], geo_col: str = "city", lat_col: str = "lat", lon_col: str = "long", n_clusters: int = 20, ) -> Dict[str, Dict[str, Any]]: """ Compare specifications with progressively richer location controls. Specifications: 1. No location controls 2. Neighbourhood / city fixed effects 3. Lat/long coordinates (linear) 4. Spatial clusters (k-means on coordinates) 5. Spatial lag of dependent variable (if possible) Parameters ---------- df : pd.DataFrame y_col : str x_cols : list[str] geo_col : str Geographic unit for fixed effects. lat_col, lon_col : str Coordinate columns. n_clusters : int Number of spatial clusters for k-means. Returns ------- dict """ import statsmodels.api as sm logger.info("Robustness: location control specifications.") sub = df.dropna(subset=[y_col] + x_cols).copy() results: Dict[str, Dict[str, Any]] = {} # 1. No location controls try: res = _safe_ols(sub[y_col], sub[x_cols]) results["No location"] = _result_to_dict(res) except Exception as e: logger.error("No-location OLS failed: %s", e) # 2. Neighbourhood / city FE if geo_col in sub.columns: try: dummies = pd.get_dummies(sub[geo_col], prefix="fe", drop_first=True, dtype=float) # Limit to top categories to avoid singular matrix top_cats = sub[geo_col].value_counts().head(30).index dummies = dummies[[c for c in dummies.columns if any(cat in c for cat in top_cats.astype(str))]] X_fe = pd.concat([sub[x_cols], dummies], axis=1) res = _safe_ols(sub[y_col], X_fe) d = _result_to_dict(res) d["FE_count"] = len(dummies.columns) results[f"{geo_col} FE"] = d except Exception as e: logger.error("FE OLS failed: %s", e) # 3. Lat/long coordinates coord_cols = [c for c in [lat_col, lon_col] if c in sub.columns] if len(coord_cols) == 2: try: sub_coord = sub.dropna(subset=coord_cols) X_coord = pd.concat( [sub_coord[x_cols], sub_coord[coord_cols]], axis=1) res = _safe_ols(sub_coord[y_col], X_coord) results["Coordinates"] = _result_to_dict(res) except Exception as e: logger.error("Coordinates OLS failed: %s", e) # 4. Spatial clusters try: from sklearn.cluster import KMeans sub_clust = sub.dropna(subset=coord_cols).copy() coords = sub_clust[coord_cols].values km = KMeans(n_clusters=min(n_clusters, len(sub_clust) // 10), random_state=42, n_init=10) sub_clust["spatial_cluster"] = km.fit_predict(coords) cl_dummies = pd.get_dummies(sub_clust["spatial_cluster"], prefix="cl", drop_first=True, dtype=float) X_cl = pd.concat([sub_clust[x_cols], cl_dummies], axis=1) res = _safe_ols(sub_clust[y_col], X_cl) results["Spatial clusters"] = _result_to_dict(res) except Exception as e: logger.error("Spatial clusters OLS failed: %s", e) # 5. Spatial lag try: from sklearn.neighbors import BallTree sub_lag = sub.dropna(subset=coord_cols).copy() tree = BallTree(np.radians(sub_lag[coord_cols].values), metric="haversine") # 5 nearest neighbours, radius in radians (~2 km) dists, indices = tree.query( np.radians(sub_lag[coord_cols].values), k=6) # Exclude self (index 0), compute mean price of neighbours neighbour_prices = [] price_arr = sub_lag[y_col].values for row_idx in range(len(indices)): nbr_idx = indices[row_idx, 1:] # skip self neighbour_prices.append(price_arr[nbr_idx].mean()) sub_lag["spatial_lag_price"] = neighbour_prices X_lag = pd.concat( [sub_lag[x_cols], sub_lag[["spatial_lag_price"]]], axis=1) res = _safe_ols(sub_lag[y_col], X_lag) results["Spatial lag"] = _result_to_dict(res) except Exception as e: logger.error("Spatial lag OLS failed: %s", e) logger.info("Location robustness: %d specs estimated.", len(results)) return results # --------------------------------------------------------------------------- # 5. Host type # --------------------------------------------------------------------------- def robustness_host_type( df: pd.DataFrame, y_col: str, x_cols: List[str], superhost_col: str = "is_superhost", listing_count_col: Optional[str] = None, min_obs: int = 50, ) -> Dict[str, Dict[str, Any]]: """ Compare regressions for different host types. Specifications: 1. All hosts 2. Single-listing hosts (if listing count available) 3. Multi-listing hosts 4. Superhosts only 5. Non-superhosts Parameters ---------- df : pd.DataFrame y_col : str x_cols : list[str] superhost_col : str listing_count_col : str, optional Column with number of listings per host. min_obs : int Returns ------- dict """ logger.info("Robustness: host type stratification.") sub = df.dropna(subset=[y_col] + x_cols).copy() results: Dict[str, Dict[str, Any]] = {} # Remove host columns from x_cols to avoid multicollinearity host_keywords = {"host", "superhost", "listing_count"} x_clean = [c for c in x_cols if not any(kw in c.lower() for kw in host_keywords)] def _run_subsample(name: str, mask: pd.Series) -> None: grp = sub.loc[mask] if len(grp) < min_obs: logger.warning("Skipping %s (n=%d).", name, len(grp)) return x_valid = [c for c in x_clean if grp[c].nunique() > 1] if not x_valid: return try: res = _safe_ols(grp[y_col], grp[x_valid]) results[name] = _result_to_dict(res) except Exception as e: logger.error("%s failed: %s", name, e) # 1. All hosts _run_subsample("All hosts", pd.Series(True, index=sub.index)) # 2/3. Single vs multi listing if listing_count_col and listing_count_col in sub.columns: _run_subsample( "Single-listing hosts", sub[listing_count_col] == 1, ) _run_subsample( "Multi-listing hosts", sub[listing_count_col] > 1, ) # 4/5. Superhosts if superhost_col in sub.columns: # Handle bool/string/numeric sh = sub[superhost_col] if sh.dtype == bool: is_super = sh elif sh.dtype == object: is_super = sh.str.lower().isin(["true", "t", "1", "yes"]) else: is_super = sh.astype(bool) _run_subsample("Superhosts", is_super) _run_subsample("Non-superhosts", ~is_super) logger.info("Host-type robustness: %d sub-samples estimated.", len(results)) return results # --------------------------------------------------------------------------- # Master runner # --------------------------------------------------------------------------- def run_all_robustness( df: pd.DataFrame, y_col: str, x_cols: List[str], room_col: str = "property_type", geo_col: str = "city", superhost_col: str = "is_superhost", listing_count_col: Optional[str] = None, ) -> Dict[str, Dict[str, Dict[str, Any]]]: """ Run all robustness checks and return a nested summary dict. Parameters ---------- df : pd.DataFrame y_col : str x_cols : list[str] room_col, geo_col, superhost_col, listing_count_col : str Returns ------- dict ``{check_category: {spec_name: results_dict}}`` """ logger.info("Running full robustness battery.") all_results: Dict[str, Dict[str, Dict[str, Any]]] = {} all_results["Winsorization"] = robustness_winsorization( df, y_col, x_cols) all_results["Room Type"] = robustness_by_room_type( df, y_col, x_cols, room_col=room_col) all_results["Price Segment"] = robustness_by_price_segment( df, y_col, x_cols) all_results["Location Controls"] = robustness_location_controls( df, y_col, x_cols, geo_col=geo_col) all_results["Host Type"] = robustness_host_type( df, y_col, x_cols, superhost_col=superhost_col, listing_count_col=listing_count_col) total_specs = sum(len(v) for v in all_results.values()) logger.info("Full robustness battery complete: %d total specifications " "across %d categories.", total_specs, len(all_results)) return all_results # --------------------------------------------------------------------------- # Summary formatter # --------------------------------------------------------------------------- def format_robustness_summary( results_dict: Dict[str, Dict[str, Dict[str, Any]]], key_variables: Optional[List[str]] = None, ) -> pd.DataFrame: """ Create a compact comparison table from all robustness results. Parameters ---------- results_dict : dict Output of ``run_all_robustness()``. key_variables : list[str], optional Subset of variables to display coefficients for. If None, only model-level statistics are shown. Returns ------- pd.DataFrame Summary table with one row per specification. """ logger.info("Formatting robustness summary.") records: List[Dict[str, Any]] = [] for category, specs in results_dict.items(): for spec_name, spec_data in specs.items(): row: Dict[str, Any] = { "Category": category, "Specification": spec_name, "N": spec_data.get("N", ""), "R2": spec_data.get("R2", ""), "Adj_R2": spec_data.get("Adj_R2", ""), "AIC": spec_data.get("AIC", ""), } # Count significant coefficients var_entries = { k: v for k, v in spec_data.items() if isinstance(v, dict) and "coef" in v } n_total = len(var_entries) n_sig = sum( 1 for v in var_entries.values() if v.get("pvalue", 1.0) < 0.05 ) row["Vars (total)"] = n_total row["Vars (p<.05)"] = n_sig # Key variable coefficients if key_variables: for kv in key_variables: if kv in spec_data and isinstance(spec_data[kv], dict): c = spec_data[kv] pval = c.get("pvalue", 1.0) stars = "" if pval <= 0.001: stars = "***" elif pval <= 0.01: stars = "**" elif pval <= 0.05: stars = "*" row[f"b({kv})"] = f"{c['coef']:.4f}{stars}" else: row[f"b({kv})"] = "" records.append(row) summary = pd.DataFrame(records) return summary