Python 64.6%
TeX 23.3%
HTML 11.6%
Makefile 0.5%
1"""Spatial econometric models for hedonic pricing analysis."""23import logging4import warnings5from typing import Dict, List, Optional, Tuple, Union67import numpy as np8import pandas as pd9import statsmodels.api as sm10from scipy.spatial import cKDTree1112logger = logging.getLogger(__name__)1314# ---------------------------------------------------------------------------15# Try importing spatial libraries (optional dependencies)16# ---------------------------------------------------------------------------1718_HAS_LIBPYSAL = False19_HAS_SPREG = False20_HAS_MGWR = False2122try:23 import libpysal24 from libpysal.weights import KNN as _KNN25 _HAS_LIBPYSAL = True26except ImportError:27 logger.info("libpysal not available; using manual KNN weights.")2829try:30 import spreg31 _HAS_SPREG = True32except ImportError:33 logger.info("spreg not available; spatial models will use 2SLS fallback.")3435try:36 import mgwr37 from mgwr.gwr import GWR as _GWR38 from mgwr.sel_bw import Sel_BW as _SelBW39 _HAS_MGWR = True40except ImportError:41 logger.info("mgwr not available; GWR will not be estimated.")424344# ---------------------------------------------------------------------------45# Helper: prepare data46# ---------------------------------------------------------------------------4748def _available_cols(df: pd.DataFrame, requested: List[str]) -> List[str]:49 """Return the subset of requested columns present in df."""50 present = [c for c in requested if c in df.columns]51 missing = set(requested) - set(present)52 if missing:53 logger.warning("Missing columns (skipped): %s", sorted(missing))54 return present555657def _prepare_data(58 df: pd.DataFrame,59 y_col: str,60 x_cols: List[str],61) -> Tuple[pd.Series, pd.DataFrame, pd.DataFrame]:62 """Drop NaN, add constant, return (y, X_with_const, clean_df)."""63 cols = _available_cols(df, x_cols)64 if not cols:65 raise ValueError("No valid feature columns found in the DataFrame.")66 lat_col = "lat" if "lat" in df.columns else "latitude"67 lon_col = "long" if "long" in df.columns else "longitude"68 need = [y_col, lat_col, lon_col] + cols69 need = list(dict.fromkeys(need)) # deduplicate, preserve order70 subset = df[need].dropna()71 if subset.empty:72 raise ValueError("No observations remaining after dropping NaN rows.")73 y = subset[y_col].astype(float)74 X = sm.add_constant(subset[cols].astype(float))75 return y, X, subset767778# ---------------------------------------------------------------------------79# Spatial weights80# ---------------------------------------------------------------------------8182class _ManualKNNWeights:83 """Minimal KNN spatial weights when libpysal is not available."""8485 def __init__(self, coords: np.ndarray, k: int = 5):86 self.n = len(coords)87 self.k = k88 tree = cKDTree(coords)89 distances, indices = tree.query(coords, k=k + 1)90 # exclude self-neighbour (index 0 in results)91 self.neighbors: Dict[int, List[int]] = {}92 self.weights: Dict[int, List[float]] = {}93 for i in range(self.n):94 nbrs = indices[i, 1:].tolist()95 self.neighbors[i] = nbrs96 self.weights[i] = [1.0 / k] * k # row-standardised9798 @property99 def sparse(self):100 """Return a scipy sparse CSR row-standardised weights matrix."""101 from scipy import sparse as sp102 rows, cols_, vals = [], [], []103 for i in range(self.n):104 for j, w in zip(self.neighbors[i], self.weights[i]):105 rows.append(i)106 cols_.append(j)107 vals.append(w)108 W = sp.csr_matrix((vals, (rows, cols_)), shape=(self.n, self.n))109 return W110111 def full(self) -> np.ndarray:112 """Dense weight matrix."""113 return self.sparse.toarray()114115116def build_weights_matrix(117 df: pd.DataFrame,118 k: int = 5,119 lat_col: str = "latitude",120 lon_col: str = "longitude",121) -> Union["_ManualKNNWeights", object]:122 """KNN-based spatial weights matrix."""123 if lat_col not in df.columns or lon_col not in df.columns:124 raise ValueError(f"DataFrame must contain '{lat_col}' and '{lon_col}' columns.")125126 coords = df[[lat_col, lon_col]].dropna().values127 n = len(coords)128 if n < k + 1:129 raise ValueError(f"Need at least k+1={k+1} observations, got {n}.")130131 if _HAS_LIBPYSAL:132 w = _KNN.from_array(coords, k=k)133 w.transform = "r" # row-standardise134 logger.info("Built libpysal KNN weights: n=%d, k=%d.", n, k)135 return w136137 w = _ManualKNNWeights(coords, k=k)138 logger.info("Built manual KNN weights: n=%d, k=%d.", n, k)139 return w140141142# ---------------------------------------------------------------------------143# Spatial autocorrelation test144# ---------------------------------------------------------------------------145146def test_spatial_autocorrelation(147 df: pd.DataFrame,148 y_col: str,149 w: object,150 permutations: int = 999,151) -> Dict[str, float]:152 """Moran's I test for spatial autocorrelation."""153 y = df[y_col].dropna().values.astype(float)154155 if _HAS_LIBPYSAL:156 from esda.moran import Moran157 mi = Moran(y, w, permutations=permutations)158 result = {159 "morans_I": float(mi.I),160 "expected_I": float(mi.EI),161 "z_score": float(mi.z_norm),162 "p_value": float(mi.p_norm),163 "p_sim": float(mi.p_sim) if hasattr(mi, "p_sim") else np.nan,164 }165 logger.info("Moran's I=%.4f, z=%.2f, p=%.4f", result["morans_I"], result["z_score"], result["p_value"])166 return result167168 # Manual Moran's I169 n = len(y)170 y_bar = y.mean()171 z = y - y_bar172173 if hasattr(w, "sparse"):174 W = w.sparse if hasattr(w.sparse, "toarray") else w.sparse175 if hasattr(W, "toarray"):176 W_dense = W.toarray()177 else:178 W_dense = np.array(W)179 elif hasattr(w, "full"):180 W_dense = w.full()181 else:182 raise ValueError("Weights object must have .sparse or .full() method.")183184 numerator = float(z @ W_dense @ z)185 denominator = float(z @ z)186 S0 = float(W_dense.sum())187188 I = (n / S0) * (numerator / denominator) if denominator != 0 else 0.0189 EI = -1.0 / (n - 1)190191 # Variance under normality assumption192 S1 = float(0.5 * ((W_dense + W_dense.T) ** 2).sum())193 S2 = float(np.sum(np.array(W_dense.sum(axis=0)).flatten() + np.array(W_dense.sum(axis=1)).flatten()) ** 2)194 # Simplified z-score195 var_I = (n ** 2 * S1 - n * S2 + 3 * S0 ** 2) / ((n ** 2 - 1) * S0 ** 2) - EI ** 2196 var_I = max(var_I, 1e-12)197 z_score = (I - EI) / np.sqrt(var_I)198199 from scipy import stats as sp_stats200 p_value = 2.0 * (1.0 - sp_stats.norm.cdf(abs(z_score)))201202 result = {203 "morans_I": float(I),204 "expected_I": float(EI),205 "z_score": float(z_score),206 "p_value": float(p_value),207 "p_sim": np.nan,208 }209 logger.info("Manual Moran's I=%.4f, z=%.2f, p=%.4f", I, z_score, p_value)210 return result211212213# ---------------------------------------------------------------------------214# Spatial lag variable creation215# ---------------------------------------------------------------------------216217def create_spatial_lags(218 df: pd.DataFrame,219 cols: List[str],220 w: object,221) -> pd.DataFrame:222 """Compute spatial lag (Wy) for specified columns."""223 result = df.copy()224 if hasattr(w, "sparse"):225 W_sp = w.sparse226 if hasattr(W_sp, "toarray"):227 W_dense = W_sp.toarray()228 else:229 W_dense = np.array(W_sp)230 elif hasattr(w, "full"):231 W_dense = w.full()232 else:233 raise ValueError("Weights object must have .sparse or .full() method.")234235 avail = _available_cols(df, cols)236 for col in avail:237 vals = df[col].values.astype(float)238 lag = W_dense @ vals239 result[f"W_{col}"] = lag240 logger.debug("Created spatial lag W_%s.", col)241242 logger.info("Created spatial lags for %d columns.", len(avail))243 return result244245246# ---------------------------------------------------------------------------247# Spatial econometric models248# ---------------------------------------------------------------------------249250def estimate_spatial_lag(251 df: pd.DataFrame,252 y_col: str,253 x_cols: List[str],254 w: object,255) -> object:256 """Spatial lag model (SAR): y = rho*Wy + Xb + e."""257 y, X, clean = _prepare_data(df, y_col, x_cols)258259 if _HAS_SPREG:260 _lat = "lat" if "lat" in clean.columns else "latitude"261 _lon = "long" if "long" in clean.columns else "longitude"262 coords = clean[[_lat, _lon]].values263 y_arr = y.values.reshape(-1, 1)264 X_arr = X.drop(columns="const", errors="ignore").values265 with warnings.catch_warnings():266 warnings.simplefilter("ignore")267 model = spreg.ML_Lag(y_arr, X_arr, w=w, name_y=y_col, name_x=list(X.columns.drop("const", errors="ignore")))268 logger.info("Spatial Lag (spreg ML): rho=%.4f, pseudo-R2=%.4f", model.rho, model.pr2)269 return model270271 # Fallback: 2SLS with Wy instrumented by WX, W2X272 return _spatial_lag_2sls(y, X, w, y_col, x_cols)273274275def _spatial_lag_2sls(276 y: pd.Series,277 X: pd.DataFrame,278 w: object,279 y_col: str,280 x_cols: List[str],281) -> object:282 """Manual 2SLS spatial lag estimation."""283 if hasattr(w, "sparse"):284 W_sp = w.sparse285 W_dense = W_sp.toarray() if hasattr(W_sp, "toarray") else np.array(W_sp)286 elif hasattr(w, "full"):287 W_dense = w.full()288 else:289 raise ValueError("Weights object must have .sparse or .full().")290291 y_arr = y.values292 Wy = W_dense @ y_arr293294 X_no_const = X.drop(columns="const", errors="ignore").values295 WX = W_dense @ X_no_const296 W2X = W_dense @ WX297298 # Instruments: X, WX, W2X299 Z = np.column_stack([np.ones(len(y_arr)), X_no_const, WX, W2X])300301 # Endogenous RHS: [const, X, Wy]302 endog_rhs = np.column_stack([np.ones(len(y_arr)), X_no_const, Wy])303304 # 2SLS305 # Stage 1: project Wy onto Z306 Pz = Z @ np.linalg.lstsq(Z, Wy, rcond=None)[0]307 # Actually need: Wy_hat = Z @ inv(Z'Z) @ Z' @ Wy308 ZtZ_inv = np.linalg.pinv(Z.T @ Z)309 Wy_hat = Z @ ZtZ_inv @ Z.T @ Wy310311 # Stage 2: regress y on [const, X, Wy_hat]312 RHS = np.column_stack([np.ones(len(y_arr)), X_no_const, Wy_hat])313 beta_2sls = np.linalg.lstsq(RHS, y_arr, rcond=None)[0]314315 # Residuals using actual Wy316 RHS_actual = np.column_stack([np.ones(len(y_arr)), X_no_const, Wy])317 resid = y_arr - RHS_actual @ beta_2sls318 n, k = RHS.shape319 sigma2 = float(resid @ resid / (n - k))320321 # Standard errors (2SLS)322 bread = np.linalg.pinv(RHS.T @ Z @ ZtZ_inv @ Z.T @ RHS)323 se = np.sqrt(np.diag(bread) * sigma2)324325 x_names = ["const"] + list(X.columns.drop("const", errors="ignore")) + ["W_" + y_col]326327 result = {328 "method": "2SLS_spatial_lag",329 "coefficients": dict(zip(x_names, beta_2sls)),330 "std_errors": dict(zip(x_names, se)),331 "rho": float(beta_2sls[-1]),332 "sigma2": sigma2,333 "n_obs": n,334 "n_params": k,335 "residuals": resid,336 "pseudo_R2": float(1.0 - np.var(resid) / np.var(y_arr)),337 }338 logger.info(339 "Spatial Lag (2SLS fallback): rho=%.4f, pseudo-R2=%.4f, n=%d",340 result["rho"], result["pseudo_R2"], n,341 )342 return result343344345def estimate_spatial_error(346 df: pd.DataFrame,347 y_col: str,348 x_cols: List[str],349 w: object,350) -> object:351 """Spatial error model (SEM): y = Xb + u, u = lambda*Wu + e."""352 y, X, clean = _prepare_data(df, y_col, x_cols)353354 if _HAS_SPREG:355 y_arr = y.values.reshape(-1, 1)356 X_arr = X.drop(columns="const", errors="ignore").values357 with warnings.catch_warnings():358 warnings.simplefilter("ignore")359 model = spreg.ML_Error(y_arr, X_arr, w=w, name_y=y_col, name_x=list(X.columns.drop("const", errors="ignore")))360 logger.info("Spatial Error (spreg ML): lambda=%.4f, pseudo-R2=%.4f", model.lam, model.pr2)361 return model362363 # Fallback: Cochrane-Orcutt-style iterative GLS364 return _spatial_error_iterative(y, X, w, y_col)365366367def _spatial_error_iterative(368 y: pd.Series,369 X: pd.DataFrame,370 w: object,371 y_col: str,372 max_iter: int = 50,373 tol: float = 1e-6,374) -> Dict:375 """Iterative estimation of spatial error model via filtered GLS."""376 if hasattr(w, "sparse"):377 W_sp = w.sparse378 W_dense = W_sp.toarray() if hasattr(W_sp, "toarray") else np.array(W_sp)379 elif hasattr(w, "full"):380 W_dense = w.full()381 else:382 raise ValueError("Weights object must have .sparse or .full().")383384 y_arr = y.values385 X_arr = X.values386 n = len(y_arr)387 I_n = np.eye(n)388389 # Initial OLS390 beta = np.linalg.lstsq(X_arr, y_arr, rcond=None)[0]391 lam = 0.0392393 for iteration in range(max_iter):394 resid = y_arr - X_arr @ beta395 Wu = W_dense @ resid396397 # Estimate lambda by regressing residuals on W*residuals398 lam_new = float(np.dot(Wu, resid) / np.dot(Wu, Wu)) if np.dot(Wu, Wu) != 0 else 0.0399 lam_new = max(-0.99, min(0.99, lam_new))400401 # Filtered variables402 A = I_n - lam_new * W_dense403 y_star = A @ y_arr404 X_star = A @ X_arr405406 beta_new = np.linalg.lstsq(X_star, y_star, rcond=None)[0]407408 if np.max(np.abs(beta_new - beta)) < tol and abs(lam_new - lam) < tol:409 beta = beta_new410 lam = lam_new411 logger.debug("Spatial error converged at iteration %d.", iteration + 1)412 break413 beta = beta_new414 lam = lam_new415416 resid_final = y_arr - X_arr @ beta417 sigma2 = float(resid_final @ resid_final / (n - len(beta)))418419 result = {420 "method": "iterative_spatial_error",421 "coefficients": dict(zip(X.columns, beta)),422 "lambda": lam,423 "sigma2": sigma2,424 "n_obs": n,425 "n_params": len(beta),426 "residuals": resid_final,427 "pseudo_R2": float(1.0 - np.var(resid_final) / np.var(y_arr)),428 }429 logger.info(430 "Spatial Error (iterative fallback): lambda=%.4f, pseudo-R2=%.4f, n=%d",431 lam, result["pseudo_R2"], n,432 )433 return result434435436def estimate_spatial_durbin(437 df: pd.DataFrame,438 y_col: str,439 x_cols: List[str],440 w: object,441) -> object:442 """Spatial Durbin model (SDM): y = rho*Wy + Xb + WX*theta + e."""443 y, X, clean = _prepare_data(df, y_col, x_cols)444445 # Create spatial lags of X446 X_no_const = X.drop(columns="const", errors="ignore")447 if hasattr(w, "sparse"):448 W_sp = w.sparse449 W_dense = W_sp.toarray() if hasattr(W_sp, "toarray") else np.array(W_sp)450 elif hasattr(w, "full"):451 W_dense = w.full()452 else:453 raise ValueError("Weights object must have .sparse or .full().")454455 WX = pd.DataFrame(456 W_dense @ X_no_const.values,457 columns=[f"W_{c}" for c in X_no_const.columns],458 index=X_no_const.index,459 )460 X_augmented = pd.concat([X_no_const, WX], axis=1)461 x_cols_augmented = list(X_augmented.columns)462463 if _HAS_SPREG:464 y_arr = y.values.reshape(-1, 1)465 X_arr = X_augmented.values466 with warnings.catch_warnings():467 warnings.simplefilter("ignore")468 model = spreg.ML_Lag(y_arr, X_arr, w=w, name_y=y_col, name_x=x_cols_augmented)469 logger.info("Spatial Durbin (spreg ML_Lag with WX): rho=%.4f, pseudo-R2=%.4f", model.rho, model.pr2)470 return model471472 # Fallback: 2SLS with augmented X473 X_aug_with_const = sm.add_constant(X_augmented)474 return _spatial_lag_2sls(y, X_aug_with_const, w, y_col, x_cols_augmented)475476477# ---------------------------------------------------------------------------478# GWR479# ---------------------------------------------------------------------------480481def run_gwr(482 df: pd.DataFrame,483 y_col: str,484 x_cols: List[str],485 lat_col: str = "latitude",486 lon_col: str = "longitude",487 kernel: str = "bisquare",488 fixed: bool = False,489) -> Optional[object]:490 """Geographically weighted regression (requires mgwr)."""491 if not _HAS_MGWR:492 logger.warning(493 "mgwr package not installed. Install with: pip install mgwr. "494 "GWR estimation is skipped."495 )496 return None497498 cols = _available_cols(df, x_cols)499 need = [y_col, lat_col, lon_col] + cols500 subset = df[need].dropna()501 if len(subset) < 50:502 logger.warning("Too few observations (%d) for reliable GWR.", len(subset))503504 coords = subset[[lat_col, lon_col]].values505 y_arr = subset[y_col].values.reshape(-1, 1)506 X_arr = subset[cols].values507508 # Bandwidth selection509 selector = _SelBW(coords, y_arr, X_arr, kernel=kernel, fixed=fixed)510 bw = selector.search()511 logger.info("GWR optimal bandwidth: %.2f", bw)512513 model = _GWR(coords, y_arr, X_arr, bw=bw, kernel=kernel, fixed=fixed)514 results = model.fit()515 logger.info(516 "GWR estimated: AICc=%.2f, R2=%.4f, n=%d",517 results.aicc, results.R2, len(y_arr),518 )519 return results520521522# ---------------------------------------------------------------------------523# Model comparison524# ---------------------------------------------------------------------------525526def compare_spatial_models(results_dict: Dict[str, object]) -> pd.DataFrame:527 """Model comparison table for spatial models."""528 rows = []529 for name, res in results_dict.items():530 row: Dict[str, object] = {"model": name}531532 if isinstance(res, dict):533 # Fallback dict results534 row["n_obs"] = res.get("n_obs", np.nan)535 row["n_params"] = res.get("n_params", np.nan)536 row["pseudo_R2"] = res.get("pseudo_R2", np.nan)537 row["sigma2"] = res.get("sigma2", np.nan)538 row["rho"] = res.get("rho", np.nan)539 row["lambda"] = res.get("lambda", np.nan)540 resid = res.get("residuals")541 if resid is not None:542 row["RMSE"] = float(np.sqrt(np.mean(resid ** 2)))543 row["MAE"] = float(np.mean(np.abs(resid)))544 row["AIC"] = _aic_from_resid(resid, res.get("n_params", 1))545 else:546 row["RMSE"] = np.nan547 row["MAE"] = np.nan548 row["AIC"] = np.nan549 elif _HAS_SPREG and hasattr(res, "pr2"):550 # spreg results551 row["n_obs"] = getattr(res, "n", np.nan)552 row["n_params"] = getattr(res, "k", np.nan)553 row["pseudo_R2"] = getattr(res, "pr2", np.nan)554 row["sigma2"] = getattr(res, "sig2", np.nan)555 row["rho"] = getattr(res, "rho", np.nan)556 row["lambda"] = getattr(res, "lam", np.nan)557 if hasattr(res, "u"):558 u = np.array(res.u).flatten()559 row["RMSE"] = float(np.sqrt(np.mean(u ** 2)))560 row["MAE"] = float(np.mean(np.abs(u)))561 row["AIC"] = _aic_from_resid(u, getattr(res, "k", 1))562 else:563 row["RMSE"] = np.nan564 row["MAE"] = np.nan565 row["AIC"] = np.nan566 else:567 # GWR or unknown568 row["n_obs"] = getattr(res, "n", np.nan)569 row["n_params"] = getattr(res, "k", np.nan)570 row["pseudo_R2"] = getattr(res, "R2", np.nan)571 row["sigma2"] = np.nan572 row["rho"] = np.nan573 row["lambda"] = np.nan574 row["RMSE"] = np.nan575 row["MAE"] = np.nan576 row["AIC"] = getattr(res, "aicc", np.nan)577578 rows.append(row)579580 comp = pd.DataFrame(rows)581 logger.info("Spatial model comparison: %d models.", len(comp))582 return comp583584585def _aic_from_resid(resid: np.ndarray, k: int) -> float:586 """Compute AIC from residuals assuming normal errors."""587 n = len(resid)588 if n <= 0:589 return np.nan590 sse = float(np.sum(resid ** 2))591 if sse <= 0:592 return np.nan593 log_lik = -0.5 * n * (np.log(2 * np.pi) + np.log(sse / n) + 1)594 return float(-2 * log_lik + 2 * k)595