Python 64.6%
TeX 23.3%
HTML 11.6%
Makefile 0.5%
1"""2robustness.py - Robustness checks for hedonic pricing analysis.34Implements a battery of specification tests following best practice in5applied micro-economics:6 - Outlier treatment (winsorisation, trimming, log, median regression)7 - Sub-sample regressions by room type and price segment8 - Progressive location control inclusion9 - Host-type stratification10 - Summary formatting for publication tables1112All functions return structured dicts that can be passed directly to13``tables.robustness_table()`` and ``tables.regression_comparison_table()``.1415Author: Simon-Pierre Boucher16Date: 2026-05-2017"""1819from __future__ import annotations2021import logging22from typing import Any, Dict, List, Optional, Tuple2324import numpy as np25import pandas as pd2627logger = logging.getLogger(__name__)282930# ---------------------------------------------------------------------------31# Internal helpers32# ---------------------------------------------------------------------------3334def _safe_ols(35 y: pd.Series,36 X: pd.DataFrame,37 add_constant: bool = True,38) -> Any:39 """Run OLS and return results; raises on failure."""40 import statsmodels.api as sm4142 if add_constant:43 X = sm.add_constant(X, has_constant="skip")44 model = sm.OLS(y, X, missing="drop")45 return model.fit(cov_type="HC1")464748def _result_to_dict(res: Any, label: str = "") -> Dict[str, Any]:49 """Convert a statsmodels result to a flat dict for comparison tables."""50 out: Dict[str, Any] = {51 "N": int(res.nobs),52 "R2": round(res.rsquared, 4),53 "Adj_R2": round(res.rsquared_adj, 4),54 "AIC": round(res.aic, 1),55 "BIC": round(res.bic, 1),56 }57 for v in res.params.index:58 out[v] = {59 "coef": round(float(res.params[v]), 6),60 "se": round(float(res.bse[v]), 6),61 "pvalue": round(float(res.pvalues[v]), 6),62 }63 return out646566def _winsorize_series(s: pd.Series, limits: Tuple[float, float] = (0.01, 0.01)) -> pd.Series:67 """Winsorize a series at the given lower/upper percentile limits."""68 from scipy.stats.mstats import winsorize as _winz69 return pd.Series(_winz(s.dropna(), limits=limits), index=s.dropna().index)707172def _trim_series(s: pd.Series, lower: float = 0.01, upper: float = 0.99) -> pd.Series:73 """Trim observations outside the specified quantile range."""74 lo = s.quantile(lower)75 hi = s.quantile(upper)76 mask = (s >= lo) & (s <= hi)77 return s[mask]787980# ---------------------------------------------------------------------------81# 1. Winsorization / outlier treatment82# ---------------------------------------------------------------------------8384def robustness_winsorization(85 df: pd.DataFrame,86 y_col: str,87 x_cols: List[str],88 limits: Tuple[float, float] = (0.01, 0.01),89 trim_quantiles: Tuple[float, float] = (0.01, 0.99),90) -> Dict[str, Dict[str, Any]]:91 """92 Compare regressions under different outlier treatments.9394 Specifications:95 1. Raw OLS96 2. Winsorized dependent variable (1st / 99th)97 3. Trimmed sample98 4. Log-transformed dependent variable99 5. Median (LAD) regression100101 Parameters102 ----------103 df : pd.DataFrame104 y_col : str105 Dependent variable (price).106 x_cols : list[str]107 Independent variables.108 limits : tuple109 Winsorization percentile limits.110 trim_quantiles : tuple111 Lower and upper quantile bounds for trimming.112113 Returns114 -------115 dict116 ``{spec_name: results_dict}``117 """118 import statsmodels.api as sm119120 logger.info("Robustness: winsorization / outlier treatment checks.")121 results: Dict[str, Dict[str, Any]] = {}122 sub = df.dropna(subset=[y_col] + x_cols).copy()123 X = sub[x_cols]124 y = sub[y_col]125126 # 1. Raw OLS127 try:128 res_raw = _safe_ols(y, X)129 results["Raw OLS"] = _result_to_dict(res_raw)130 except Exception as e:131 logger.error("Raw OLS failed: %s", e)132133 # 2. Winsorized134 try:135 y_wins = _winsorize_series(y, limits=limits)136 # Align X with winsorized y index137 X_wins = X.loc[y_wins.index]138 res_wins = _safe_ols(y_wins, X_wins)139 results["Winsorized (1%/99%)"] = _result_to_dict(res_wins)140 except Exception as e:141 logger.error("Winsorized OLS failed: %s", e)142143 # 3. Trimmed144 try:145 lo = y.quantile(trim_quantiles[0])146 hi = y.quantile(trim_quantiles[1])147 mask = (y >= lo) & (y <= hi)148 res_trim = _safe_ols(y[mask], X.loc[mask])149 results["Trimmed (1%-99%)"] = _result_to_dict(res_trim)150 except Exception as e:151 logger.error("Trimmed OLS failed: %s", e)152153 # 4. Log-transformed154 try:155 pos_mask = y > 0156 y_log = np.log(y[pos_mask])157 res_log = _safe_ols(y_log, X.loc[pos_mask])158 results["Log(price) OLS"] = _result_to_dict(res_log)159 except Exception as e:160 logger.error("Log OLS failed: %s", e)161162 # 5. Median regression (LAD)163 try:164 X_const = sm.add_constant(X, has_constant="skip")165 model_med = sm.QuantReg(y, X_const)166 res_med = model_med.fit(q=0.5)167 med_dict: Dict[str, Any] = {168 "N": int(res_med.nobs),169 "R2": np.nan, # pseudo R2 not directly comparable170 "Adj_R2": np.nan,171 }172 # Pseudo R2 for quantile regression173 try:174 med_dict["Pseudo_R2"] = round(float(1 - res_med.prsquared), 4)175 except Exception:176 pass177 for v in res_med.params.index:178 med_dict[v] = {179 "coef": round(float(res_med.params[v]), 6),180 "se": round(float(res_med.bse[v]), 6),181 "pvalue": round(float(res_med.pvalues[v]), 6),182 }183 results["Median Regression"] = med_dict184 except Exception as e:185 logger.error("Median regression failed: %s", e)186187 logger.info("Winsorization robustness: %d/%d specs succeeded.",188 len(results), 5)189 return results190191192# ---------------------------------------------------------------------------193# 2. By room type194# ---------------------------------------------------------------------------195196def robustness_by_room_type(197 df: pd.DataFrame,198 y_col: str,199 x_cols: List[str],200 room_col: str = "property_type",201 min_obs: int = 50,202) -> Dict[str, Dict[str, Any]]:203 """204 Separate regressions for each room/property type.205206 Parameters207 ----------208 df : pd.DataFrame209 y_col : str210 x_cols : list[str]211 room_col : str212 Column identifying room/property type.213 min_obs : int214 Minimum observations to run a regression for a given type.215216 Returns217 -------218 dict219 ``{room_type: results_dict}``220 """221 logger.info("Robustness: sub-sample by room type (%s).", room_col)222 sub = df.dropna(subset=[y_col] + x_cols + [room_col]).copy()223 results: Dict[str, Dict[str, Any]] = {}224225 # Full sample first226 try:227 res_full = _safe_ols(sub[y_col], sub[x_cols])228 results["Full Sample"] = _result_to_dict(res_full)229 except Exception as e:230 logger.error("Full sample OLS failed: %s", e)231232 for rtype, grp in sub.groupby(room_col, observed=True):233 if len(grp) < min_obs:234 logger.warning("Skipping %s (n=%d < %d).", rtype, len(grp), min_obs)235 continue236 # Remove x_cols that have zero variance in this sub-sample237 x_valid = [c for c in x_cols if grp[c].nunique() > 1]238 if not x_valid:239 logger.warning("Skipping %s: no valid regressors.", rtype)240 continue241 try:242 res = _safe_ols(grp[y_col], grp[x_valid])243 results[str(rtype)] = _result_to_dict(res)244 except Exception as e:245 logger.error("OLS for %s failed: %s", rtype, e)246247 logger.info("Room-type robustness: %d sub-samples estimated.", len(results))248 return results249250251# ---------------------------------------------------------------------------252# 3. By price segment253# ---------------------------------------------------------------------------254255def robustness_by_price_segment(256 df: pd.DataFrame,257 y_col: str,258 x_cols: List[str],259 segments: Optional[Dict[str, Tuple[float, float]]] = None,260 price_category_col: Optional[str] = "price_category",261 min_obs: int = 50,262) -> Dict[str, Dict[str, Any]]:263 """264 Separate regressions for budget, moderate, premium, luxury segments.265266 Parameters267 ----------268 df : pd.DataFrame269 y_col : str270 x_cols : list[str]271 segments : dict, optional272 ``{segment_name: (lower_bound, upper_bound)}``.273 If None, quartile-based segments are used.274 price_category_col : str, optional275 If present in df, use this column instead of numeric cuts.276 min_obs : int277 Minimum observations per segment.278279 Returns280 -------281 dict282 """283 logger.info("Robustness: sub-sample by price segment.")284 sub = df.dropna(subset=[y_col] + x_cols).copy()285 results: Dict[str, Dict[str, Any]] = {}286287 # Full sample288 try:289 res_full = _safe_ols(sub[y_col], sub[x_cols])290 results["Full Sample"] = _result_to_dict(res_full)291 except Exception as e:292 logger.error("Full sample OLS failed: %s", e)293294 if price_category_col and price_category_col in sub.columns:295 for seg, grp in sub.groupby(price_category_col, observed=True):296 if len(grp) < min_obs:297 continue298 x_valid = [c for c in x_cols if grp[c].nunique() > 1]299 if not x_valid:300 continue301 try:302 res = _safe_ols(grp[y_col], grp[x_valid])303 results[str(seg)] = _result_to_dict(res)304 except Exception as e:305 logger.error("Segment %s failed: %s", seg, e)306 else:307 if segments is None:308 q25, q50, q75 = sub[y_col].quantile([0.25, 0.50, 0.75])309 segments = {310 "Budget (Q1)": (sub[y_col].min(), q25),311 "Moderate (Q2)": (q25, q50),312 "Premium (Q3)": (q50, q75),313 "Luxury (Q4)": (q75, sub[y_col].max()),314 }315316 for seg_name, (lo, hi) in segments.items():317 grp = sub[(sub[y_col] >= lo) & (sub[y_col] <= hi)]318 if len(grp) < min_obs:319 logger.warning("Skipping %s (n=%d).", seg_name, len(grp))320 continue321 x_valid = [c for c in x_cols if grp[c].nunique() > 1]322 if not x_valid:323 continue324 try:325 res = _safe_ols(grp[y_col], grp[x_valid])326 results[seg_name] = _result_to_dict(res)327 except Exception as e:328 logger.error("Segment %s failed: %s", seg_name, e)329330 logger.info("Price-segment robustness: %d segments estimated.", len(results))331 return results332333334# ---------------------------------------------------------------------------335# 4. Location controls336# ---------------------------------------------------------------------------337338def robustness_location_controls(339 df: pd.DataFrame,340 y_col: str,341 x_cols: List[str],342 geo_col: str = "city",343 lat_col: str = "lat",344 lon_col: str = "long",345 n_clusters: int = 20,346) -> Dict[str, Dict[str, Any]]:347 """348 Compare specifications with progressively richer location controls.349350 Specifications:351 1. No location controls352 2. Neighbourhood / city fixed effects353 3. Lat/long coordinates (linear)354 4. Spatial clusters (k-means on coordinates)355 5. Spatial lag of dependent variable (if possible)356357 Parameters358 ----------359 df : pd.DataFrame360 y_col : str361 x_cols : list[str]362 geo_col : str363 Geographic unit for fixed effects.364 lat_col, lon_col : str365 Coordinate columns.366 n_clusters : int367 Number of spatial clusters for k-means.368369 Returns370 -------371 dict372 """373 import statsmodels.api as sm374375 logger.info("Robustness: location control specifications.")376 sub = df.dropna(subset=[y_col] + x_cols).copy()377 results: Dict[str, Dict[str, Any]] = {}378379 # 1. No location controls380 try:381 res = _safe_ols(sub[y_col], sub[x_cols])382 results["No location"] = _result_to_dict(res)383 except Exception as e:384 logger.error("No-location OLS failed: %s", e)385386 # 2. Neighbourhood / city FE387 if geo_col in sub.columns:388 try:389 dummies = pd.get_dummies(sub[geo_col], prefix="fe", drop_first=True,390 dtype=float)391 # Limit to top categories to avoid singular matrix392 top_cats = sub[geo_col].value_counts().head(30).index393 dummies = dummies[[c for c in dummies.columns394 if any(cat in c for cat in top_cats.astype(str))]]395 X_fe = pd.concat([sub[x_cols], dummies], axis=1)396 res = _safe_ols(sub[y_col], X_fe)397 d = _result_to_dict(res)398 d["FE_count"] = len(dummies.columns)399 results[f"{geo_col} FE"] = d400 except Exception as e:401 logger.error("FE OLS failed: %s", e)402403 # 3. Lat/long coordinates404 coord_cols = [c for c in [lat_col, lon_col] if c in sub.columns]405 if len(coord_cols) == 2:406 try:407 sub_coord = sub.dropna(subset=coord_cols)408 X_coord = pd.concat(409 [sub_coord[x_cols], sub_coord[coord_cols]], axis=1)410 res = _safe_ols(sub_coord[y_col], X_coord)411 results["Coordinates"] = _result_to_dict(res)412 except Exception as e:413 logger.error("Coordinates OLS failed: %s", e)414415 # 4. Spatial clusters416 try:417 from sklearn.cluster import KMeans418 sub_clust = sub.dropna(subset=coord_cols).copy()419 coords = sub_clust[coord_cols].values420 km = KMeans(n_clusters=min(n_clusters, len(sub_clust) // 10),421 random_state=42, n_init=10)422 sub_clust["spatial_cluster"] = km.fit_predict(coords)423 cl_dummies = pd.get_dummies(sub_clust["spatial_cluster"],424 prefix="cl", drop_first=True,425 dtype=float)426 X_cl = pd.concat([sub_clust[x_cols], cl_dummies], axis=1)427 res = _safe_ols(sub_clust[y_col], X_cl)428 results["Spatial clusters"] = _result_to_dict(res)429 except Exception as e:430 logger.error("Spatial clusters OLS failed: %s", e)431432 # 5. Spatial lag433 try:434 from sklearn.neighbors import BallTree435 sub_lag = sub.dropna(subset=coord_cols).copy()436 tree = BallTree(np.radians(sub_lag[coord_cols].values),437 metric="haversine")438 # 5 nearest neighbours, radius in radians (~2 km)439 dists, indices = tree.query(440 np.radians(sub_lag[coord_cols].values), k=6)441 # Exclude self (index 0), compute mean price of neighbours442 neighbour_prices = []443 price_arr = sub_lag[y_col].values444 for row_idx in range(len(indices)):445 nbr_idx = indices[row_idx, 1:] # skip self446 neighbour_prices.append(price_arr[nbr_idx].mean())447 sub_lag["spatial_lag_price"] = neighbour_prices448 X_lag = pd.concat(449 [sub_lag[x_cols],450 sub_lag[["spatial_lag_price"]]], axis=1)451 res = _safe_ols(sub_lag[y_col], X_lag)452 results["Spatial lag"] = _result_to_dict(res)453 except Exception as e:454 logger.error("Spatial lag OLS failed: %s", e)455456 logger.info("Location robustness: %d specs estimated.", len(results))457 return results458459460# ---------------------------------------------------------------------------461# 5. Host type462# ---------------------------------------------------------------------------463464def robustness_host_type(465 df: pd.DataFrame,466 y_col: str,467 x_cols: List[str],468 superhost_col: str = "is_superhost",469 listing_count_col: Optional[str] = None,470 min_obs: int = 50,471) -> Dict[str, Dict[str, Any]]:472 """473 Compare regressions for different host types.474475 Specifications:476 1. All hosts477 2. Single-listing hosts (if listing count available)478 3. Multi-listing hosts479 4. Superhosts only480 5. Non-superhosts481482 Parameters483 ----------484 df : pd.DataFrame485 y_col : str486 x_cols : list[str]487 superhost_col : str488 listing_count_col : str, optional489 Column with number of listings per host.490 min_obs : int491492 Returns493 -------494 dict495 """496 logger.info("Robustness: host type stratification.")497 sub = df.dropna(subset=[y_col] + x_cols).copy()498 results: Dict[str, Dict[str, Any]] = {}499500 # Remove host columns from x_cols to avoid multicollinearity501 host_keywords = {"host", "superhost", "listing_count"}502 x_clean = [c for c in x_cols503 if not any(kw in c.lower() for kw in host_keywords)]504505 def _run_subsample(name: str, mask: pd.Series) -> None:506 grp = sub.loc[mask]507 if len(grp) < min_obs:508 logger.warning("Skipping %s (n=%d).", name, len(grp))509 return510 x_valid = [c for c in x_clean if grp[c].nunique() > 1]511 if not x_valid:512 return513 try:514 res = _safe_ols(grp[y_col], grp[x_valid])515 results[name] = _result_to_dict(res)516 except Exception as e:517 logger.error("%s failed: %s", name, e)518519 # 1. All hosts520 _run_subsample("All hosts", pd.Series(True, index=sub.index))521522 # 2/3. Single vs multi listing523 if listing_count_col and listing_count_col in sub.columns:524 _run_subsample(525 "Single-listing hosts",526 sub[listing_count_col] == 1,527 )528 _run_subsample(529 "Multi-listing hosts",530 sub[listing_count_col] > 1,531 )532533 # 4/5. Superhosts534 if superhost_col in sub.columns:535 # Handle bool/string/numeric536 sh = sub[superhost_col]537 if sh.dtype == bool:538 is_super = sh539 elif sh.dtype == object:540 is_super = sh.str.lower().isin(["true", "t", "1", "yes"])541 else:542 is_super = sh.astype(bool)543544 _run_subsample("Superhosts", is_super)545 _run_subsample("Non-superhosts", ~is_super)546547 logger.info("Host-type robustness: %d sub-samples estimated.", len(results))548 return results549550551# ---------------------------------------------------------------------------552# Master runner553# ---------------------------------------------------------------------------554555def run_all_robustness(556 df: pd.DataFrame,557 y_col: str,558 x_cols: List[str],559 room_col: str = "property_type",560 geo_col: str = "city",561 superhost_col: str = "is_superhost",562 listing_count_col: Optional[str] = None,563) -> Dict[str, Dict[str, Dict[str, Any]]]:564 """565 Run all robustness checks and return a nested summary dict.566567 Parameters568 ----------569 df : pd.DataFrame570 y_col : str571 x_cols : list[str]572 room_col, geo_col, superhost_col, listing_count_col : str573574 Returns575 -------576 dict577 ``{check_category: {spec_name: results_dict}}``578 """579 logger.info("Running full robustness battery.")580 all_results: Dict[str, Dict[str, Dict[str, Any]]] = {}581582 all_results["Winsorization"] = robustness_winsorization(583 df, y_col, x_cols)584585 all_results["Room Type"] = robustness_by_room_type(586 df, y_col, x_cols, room_col=room_col)587588 all_results["Price Segment"] = robustness_by_price_segment(589 df, y_col, x_cols)590591 all_results["Location Controls"] = robustness_location_controls(592 df, y_col, x_cols, geo_col=geo_col)593594 all_results["Host Type"] = robustness_host_type(595 df, y_col, x_cols, superhost_col=superhost_col,596 listing_count_col=listing_count_col)597598 total_specs = sum(len(v) for v in all_results.values())599 logger.info("Full robustness battery complete: %d total specifications "600 "across %d categories.", total_specs, len(all_results))601 return all_results602603604# ---------------------------------------------------------------------------605# Summary formatter606# ---------------------------------------------------------------------------607608def format_robustness_summary(609 results_dict: Dict[str, Dict[str, Dict[str, Any]]],610 key_variables: Optional[List[str]] = None,611) -> pd.DataFrame:612 """613 Create a compact comparison table from all robustness results.614615 Parameters616 ----------617 results_dict : dict618 Output of ``run_all_robustness()``.619 key_variables : list[str], optional620 Subset of variables to display coefficients for.621 If None, only model-level statistics are shown.622623 Returns624 -------625 pd.DataFrame626 Summary table with one row per specification.627 """628 logger.info("Formatting robustness summary.")629 records: List[Dict[str, Any]] = []630631 for category, specs in results_dict.items():632 for spec_name, spec_data in specs.items():633 row: Dict[str, Any] = {634 "Category": category,635 "Specification": spec_name,636 "N": spec_data.get("N", ""),637 "R2": spec_data.get("R2", ""),638 "Adj_R2": spec_data.get("Adj_R2", ""),639 "AIC": spec_data.get("AIC", ""),640 }641642 # Count significant coefficients643 var_entries = {644 k: v for k, v in spec_data.items()645 if isinstance(v, dict) and "coef" in v646 }647 n_total = len(var_entries)648 n_sig = sum(649 1 for v in var_entries.values()650 if v.get("pvalue", 1.0) < 0.05651 )652 row["Vars (total)"] = n_total653 row["Vars (p<.05)"] = n_sig654655 # Key variable coefficients656 if key_variables:657 for kv in key_variables:658 if kv in spec_data and isinstance(spec_data[kv], dict):659 c = spec_data[kv]660 pval = c.get("pvalue", 1.0)661 stars = ""662 if pval <= 0.001:663 stars = "***"664 elif pval <= 0.01:665 stars = "**"666 elif pval <= 0.05:667 stars = "*"668 row[f"b({kv})"] = f"{c['coef']:.4f}{stars}"669 else:670 row[f"b({kv})"] = ""671672 records.append(row)673674 summary = pd.DataFrame(records)675 return summary676