Python 64.6%
TeX 23.3%
HTML 11.6%
Makefile 0.5%
1"""Core hedonic regression models for Airbnb pricing analysis."""23import logging4from typing import Dict, List, Optional, Tuple56import numpy as np7import pandas as pd8import statsmodels.api as sm9from sklearn.model_selection import KFold10from sklearn.metrics import mean_squared_error, mean_absolute_error1112logger = logging.getLogger(__name__)1314# ---------------------------------------------------------------------------15# Feature group definitions for the 9-model nested sequence16# ---------------------------------------------------------------------------1718STRUCTURAL_COLS = [19 "bedrooms", "bathrooms", "guests_count", "is_entire_home",20 "is_private_room", "log_accommodates", "log_bedrooms", "log_bathrooms",21 "bathrooms_per_guest", "bedrooms_per_guest", "large_capacity_listing",22 "accommodates_sq", "bedrooms_sq", "bathrooms_sq",23]2425AMENITY_COLS = [26 "amenities_count", "amenity_richness_index", "luxury_index",27 "business_travel_index", "family_comfort_index", "leisure_index",28 "remote_work_index", "pets_allowed",29]3031HOST_COLS = [32 "is_superhost", "superhost_flag", "host_professional_dummy",33 "multi_listing_host", "superhost_x_entire_home", "guest_favorite_flag",34 "is_guest_favorite",35]3637REVIEW_COLS = [38 "rating", "num_reviews", "log_reviews", "review_quality",39 "reputation_score", "cleanliness_premium", "location_rating_premium",40 "has_reviews", "quality_score",41]4243LOCATION_COLS = [44 "distance_to_center", "lat", "long",45]4647NEIGHBOURHOOD_COLS = [48 "spatial_cluster",49]5051NLP_COLS = [52 "title_length", "title_word_count", "luxury_language_score",53 "cozy_language_score", "professional_tone_score", "tourism_language_score",54 "text_embed_0", "text_embed_1", "text_embed_2", "text_embed_3",55]5657SPATIAL_CONTROL_COLS = [58 "local_density_500m", "local_density_1km",59 "local_avg_price_500m", "local_avg_price_1km",60]6162MODEL_SEQUENCE = {63 "Model 1 - Structural": STRUCTURAL_COLS,64 "Model 2 - + Amenities": STRUCTURAL_COLS + AMENITY_COLS,65 "Model 3 - + Host": STRUCTURAL_COLS + AMENITY_COLS + HOST_COLS,66 "Model 4 - + Reviews": STRUCTURAL_COLS + AMENITY_COLS + HOST_COLS + REVIEW_COLS,67 "Model 5 - + Location": STRUCTURAL_COLS + AMENITY_COLS + HOST_COLS + REVIEW_COLS + LOCATION_COLS,68 "Model 6 - + Neighbourhood": STRUCTURAL_COLS + AMENITY_COLS + HOST_COLS + REVIEW_COLS + LOCATION_COLS + NEIGHBOURHOOD_COLS,69 "Model 7 - + NLP": STRUCTURAL_COLS + AMENITY_COLS + HOST_COLS + REVIEW_COLS + LOCATION_COLS + NEIGHBOURHOOD_COLS + NLP_COLS,70 "Model 8 - + Spatial": STRUCTURAL_COLS + AMENITY_COLS + HOST_COLS + REVIEW_COLS + LOCATION_COLS + NEIGHBOURHOOD_COLS + NLP_COLS + SPATIAL_CONTROL_COLS,71 "Model 9 - Full Hedonic": STRUCTURAL_COLS + AMENITY_COLS + HOST_COLS + REVIEW_COLS + LOCATION_COLS + NEIGHBOURHOOD_COLS + NLP_COLS + SPATIAL_CONTROL_COLS,72}737475def _available_cols(df: pd.DataFrame, requested: List[str]) -> List[str]:76 """Return the subset of requested columns that actually exist in df."""77 present = [c for c in requested if c in df.columns]78 missing = set(requested) - set(present)79 if missing:80 logger.warning("Missing columns (will be skipped): %s", sorted(missing))81 return present828384def _prepare_data(85 df: pd.DataFrame,86 y_col: str,87 x_cols: List[str],88) -> Tuple[pd.Series, pd.DataFrame]:89 """Drop rows with NaN in y or X; add constant; return (y, X)."""90 cols = _available_cols(df, x_cols)91 if not cols:92 raise ValueError("No valid feature columns found in the DataFrame.")93 subset = df[[y_col] + cols].dropna()94 if subset.empty:95 raise ValueError("No observations remaining after dropping NaN rows.")96 y = subset[y_col].astype(float)97 X = sm.add_constant(subset[cols].astype(float))98 logger.info("Prepared data: n=%d, k=%d (excl. constant)", len(y), len(cols))99 return y, X100101102# ---------------------------------------------------------------------------103# OLS estimators104# ---------------------------------------------------------------------------105106def estimate_ols(107 df: pd.DataFrame,108 y_col: str,109 x_cols: List[str],110 robust: str = "HC3",111) -> sm.regression.linear_model.RegressionResultsWrapper:112 """OLS with heteroskedasticity-robust standard errors."""113 y, X = _prepare_data(df, y_col, x_cols)114 model = sm.OLS(y, X)115 results = model.fit(cov_type=robust)116 logger.info(117 "OLS (%s): R2=%.4f, Adj-R2=%.4f, n=%d",118 robust, results.rsquared, results.rsquared_adj, int(results.nobs),119 )120 return results121122123def estimate_ols_clustered(124 df: pd.DataFrame,125 y_col: str,126 x_cols: List[str],127 cluster_col: str,128) -> sm.regression.linear_model.RegressionResultsWrapper:129 """OLS with cluster-robust standard errors."""130 cols = _available_cols(df, x_cols)131 subset = df[[y_col] + cols + [cluster_col]].dropna()132 if subset.empty:133 raise ValueError("No observations remaining after dropping NaN rows.")134 y = subset[y_col].astype(float)135 X = sm.add_constant(subset[cols].astype(float))136 groups = subset[cluster_col]137 model = sm.OLS(y, X)138 results = model.fit(cov_type="cluster", cov_kwds={"groups": groups})139 logger.info(140 "OLS (clustered on %s): R2=%.4f, n=%d, clusters=%d",141 cluster_col, results.rsquared, int(results.nobs), groups.nunique(),142 )143 return results144145146# ---------------------------------------------------------------------------147# Nested model sequence148# ---------------------------------------------------------------------------149150def build_model_sequence(151 df: pd.DataFrame,152 y_col: str = "log_price",153 robust: str = "HC3",154) -> Dict[str, sm.regression.linear_model.RegressionResultsWrapper]:155 """Estimate the 9-model nested hedonic sequence."""156 results: Dict[str, sm.regression.linear_model.RegressionResultsWrapper] = {}157 for name, xcols in MODEL_SEQUENCE.items():158 avail = _available_cols(df, xcols)159 if not avail:160 logger.warning("Skipping %s: no available features.", name)161 continue162 try:163 res = estimate_ols(df, y_col, avail, robust=robust)164 results[name] = res165 logger.info("Estimated %s (k=%d)", name, len(avail))166 except Exception as exc:167 logger.error("Failed to estimate %s: %s", name, exc)168 return results169170171# ---------------------------------------------------------------------------172# Model comparison & diagnostics173# ---------------------------------------------------------------------------174175def compare_models(176 results_dict: Dict[str, sm.regression.linear_model.RegressionResultsWrapper],177) -> pd.DataFrame:178 """Extract R2, adj-R2, RMSE, MAE, AIC, BIC for each model."""179 rows = []180 for name, res in results_dict.items():181 resid = res.resid182 n = int(res.nobs)183 rmse = np.sqrt(np.mean(resid ** 2))184 mae = np.mean(np.abs(resid))185 rows.append({186 "model": name,187 "n_obs": n,188 "n_params": int(res.df_model) + 1,189 "R2": round(res.rsquared, 6),190 "Adj_R2": round(res.rsquared_adj, 6),191 "RMSE": round(rmse, 6),192 "MAE": round(mae, 6),193 "AIC": round(res.aic, 2),194 "BIC": round(res.bic, 2),195 })196 comp = pd.DataFrame(rows)197 logger.info("Model comparison table built (%d models).", len(comp))198 return comp199200201def compute_price_premiums(202 results: sm.regression.linear_model.RegressionResultsWrapper,203) -> pd.DataFrame:204 """Convert log-linear coefficients to percentage price premiums: 100*(exp(b)-1)."""205 params = results.params.drop("const", errors="ignore")206 se = results.bse.drop("const", errors="ignore")207 pvals = results.pvalues.drop("const", errors="ignore")208 ci = results.conf_int().drop("const", errors="ignore")209210 premium = 100.0 * (np.exp(params) - 1)211 premium_lo = 100.0 * (np.exp(ci[0]) - 1)212 premium_hi = 100.0 * (np.exp(ci[1]) - 1)213214 df_out = pd.DataFrame({215 "coefficient": params.values,216 "std_error": se.values,217 "p_value": pvals.values,218 "premium_pct": premium.values,219 "premium_lo_95": premium_lo.values,220 "premium_hi_95": premium_hi.values,221 }, index=params.index)222 df_out = df_out.sort_values("premium_pct", ascending=False)223 logger.info("Computed price premiums for %d variables.", len(df_out))224 return df_out225226227# ---------------------------------------------------------------------------228# Cross-validation229# ---------------------------------------------------------------------------230231def cross_validate_ols(232 df: pd.DataFrame,233 y_col: str,234 x_cols: List[str],235 k: int = 5,236 random_state: int = 42,237) -> Dict[str, float]:238 """K-fold cross-validation returning out-of-sample RMSE and MAE."""239 y, X = _prepare_data(df, y_col, x_cols)240 kf = KFold(n_splits=k, shuffle=True, random_state=random_state)241 rmses, maes = [], []242 for fold, (train_idx, test_idx) in enumerate(kf.split(X), 1):243 X_train, X_test = X.iloc[train_idx], X.iloc[test_idx]244 y_train, y_test = y.iloc[train_idx], y.iloc[test_idx]245 model = sm.OLS(y_train, X_train).fit()246 y_pred = model.predict(X_test)247 rmses.append(np.sqrt(mean_squared_error(y_test, y_pred)))248 maes.append(mean_absolute_error(y_test, y_pred))249 logger.debug("Fold %d: RMSE=%.4f, MAE=%.4f", fold, rmses[-1], maes[-1])250251 result = {252 "mean_rmse": float(np.mean(rmses)),253 "std_rmse": float(np.std(rmses)),254 "mean_mae": float(np.mean(maes)),255 "std_mae": float(np.std(maes)),256 "k": k,257 }258 logger.info(259 "%d-fold CV: RMSE=%.4f (+/-%.4f), MAE=%.4f (+/-%.4f)",260 k, result["mean_rmse"], result["std_rmse"],261 result["mean_mae"], result["std_mae"],262 )263 return result264265266# ---------------------------------------------------------------------------267# Publication-quality regression table268# ---------------------------------------------------------------------------269270def format_regression_table(271 results_dict: Dict[str, sm.regression.linear_model.RegressionResultsWrapper],272 stars: bool = True,273) -> pd.DataFrame:274 """Create a publication-ready regression comparison table (coef + SE rows)."""275 all_vars: List[str] = []276 for res in results_dict.values():277 for v in res.params.index:278 if v not in all_vars:279 all_vars.append(v)280281 rows: List[Dict[str, str]] = []282 for var in all_vars:283 coef_row: Dict[str, str] = {"variable": var}284 se_row: Dict[str, str] = {"variable": ""}285 for mname, res in results_dict.items():286 if var in res.params.index:287 b = res.params[var]288 s = res.bse[var]289 p = res.pvalues[var]290 star_str = ""291 if stars:292 if p < 0.001:293 star_str = "***"294 elif p < 0.01:295 star_str = "**"296 elif p < 0.05:297 star_str = "*"298 coef_row[mname] = f"{b:.4f}{star_str}"299 se_row[mname] = f"({s:.4f})"300 else:301 coef_row[mname] = ""302 se_row[mname] = ""303 rows.append(coef_row)304 rows.append(se_row)305306 # Footer rows307 for stat_name, extractor in [308 ("N", lambda r: f"{int(r.nobs)}"),309 ("R2", lambda r: f"{r.rsquared:.4f}"),310 ("Adj. R2", lambda r: f"{r.rsquared_adj:.4f}"),311 ("AIC", lambda r: f"{r.aic:.1f}"),312 ("BIC", lambda r: f"{r.bic:.1f}"),313 ]:314 row: Dict[str, str] = {"variable": stat_name}315 for mname, res in results_dict.items():316 row[mname] = extractor(res)317 rows.append(row)318319 table = pd.DataFrame(rows)320 logger.info("Formatted regression table: %d rows x %d models.", len(table), len(results_dict))321 return table322