xjtulyc Housing Market Analysis
When to Use This Skill
Use this skill when you need to:
- Estimate implicit prices of housing attributes (hedonic regression)
- Construct repeat-sales house price indices (Case-Shiller method)
- Measure housing affordability and price-to-income ratios
- Detect spatial clustering of housing prices (Moran's I, LISA)
- Forecast house prices using machine learning and spatial lag models
- Analyze rental yields and capitalization rates
- Study filtering theory and housing supply/demand dynamics
Trigger keywords: hedonic pricing, house price index, repeat sales, affordability ratio, housing bubble, capitalization rate, rental yield, spatial lag model, housing supply, zoning analysis, price-to-rent, mortgage payment, housing filtering, location gradient, CBD distance gradient, amenity valuation.
Background & Key Concepts
Hedonic Price Model
Houses are bundles of attributes. The hedonic regression decomposes price $P$ into contributions of individual characteristics:
$$\ln P_i = \alpha + \sum_k \beta_k X_{ik} + \varepsilon_i$$
where $X_{ik}$ includes structural attributes (sqft, bedrooms, age), location variables (distance to CBD, school quality), and neighborhood characteristics. The implicit price of attribute $k$ is $\partial P/\partial X_k = \beta_k \cdot P$.
Repeat-Sales Index
For two sales of property $i$ at times $s$ and $t$:
$$\ln P_{it} - \ln P_{is} = \sum_{\tau=s+1}^{t} \delta_\tau D_\tau + \varepsilon_{it}$$
where $D_\tau$ are time dummies and $\delta_\tau$ estimates log price changes. The Case-Shiller weighted repeat-sales accounts for variance growing with holding period.
Affordability Measures
- Price-to-Income Ratio (PIR): $\text{PIR} = P / Y_{\text{median}}$
- Housing Affordability Index (HAI): $\text{HAI} = \frac{\text{Qualifying Income}}{Y_{\text{median}}} \times 100$, where qualifying income = monthly payment × 12 / 0.28
- Residual Income: $I - \text{Housing Cost} - \text{Transport Cost}$ (location efficiency)
Spatial Hedonic Model
$$\ln \mathbf{P} = \rho \mathbf{W} \ln \mathbf{P} + \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$$
The spatial autoregressive (SAR) model accounts for spatial spillovers: $\rho$ captures the degree to which nearby prices influence each other.
Von Thünen / Bid-Rent Gradient
$$P(d) = P_0 \exp(-\gamma d)$$
Semilog regression of $\ln P$ on distance $d$ from CBD gives gradient $\gamma$ (price decay per km).
Environment Setup
pip install numpy>=1.24 pandas>=2.0 scipy>=1.11 statsmodels>=0.14 \
scikit-learn>=1.3 matplotlib>=3.7 geopandas>=0.14
import numpy as np
import pandas as pd
import statsmodels.api as sm
from sklearn.ensemble import GradientBoostingRegressor
import matplotlib.pyplot as plt
print("Housing market analysis environment ready")
Core Workflow
Step 1: Hedonic Regression with Spatial Effects
import numpy as np
import pandas as pd
import statsmodels.api as sm
from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt
# -----------------------------------------------------------------
# Simulate hedonic dataset: 1000 transactions
# -----------------------------------------------------------------
np.random.seed(42)
n = 1000
# Location: x, y coordinates in a 20km x 20km city
x = np.random.uniform(0, 20, n)
y = np.random.uniform(0, 20, n)
# Distance to CBD at (10, 10)
cbd_dist = np.sqrt((x - 10)**2 + (y - 10)**2)
# Structural attributes
sqft = np.random.lognormal(7.0, 0.4, n) # 330-2700 sqft
bedrooms = np.random.choice([1, 2, 3, 4, 5], n, p=[0.1, 0.25, 0.35, 0.2, 0.1])
age = np.random.randint(0, 80, n).astype(float)
garage = np.random.binomial(1, 0.6, n)
school_rating = np.random.uniform(3, 10, n)
# True hedonic model (log-linear)
ln_P = (
11.0 # intercept ~ $60k baseline
+ 0.7 * np.log(sqft) # size elasticity 0.7
+ 0.05 * bedrooms # each bedroom +5%
- 0.005 * age # depreciation
+ 0.10 * garage # garage premium
+ 0.06 * school_rating # school quality
- 0.04 * cbd_dist # CBD gradient: -4%/km
+ np.random.normal(0, 0.12, n) # unexplained
)
price = np.exp(ln_P)
df = pd.DataFrame({
"price": price, "sqft": sqft, "bedrooms": bedrooms,
"age": age, "garage": garage, "school_rating": school_rating,
"cbd_dist": cbd_dist, "x": x, "y": y
})
# -----------------------------------------------------------------
# OLS hedonic regression
# -----------------------------------------------------------------
df["ln_price"] = np.log(df["price"])
df["ln_sqft"] = np.log(df["sqft"])
features = ["ln_sqft", "bedrooms", "age", "garage", "school_rating", "cbd_dist"]
X_reg = sm.add_constant(df[features])
ols = sm.OLS(df["ln_price"], X_reg).fit(cov_type="HC3")
print("=== OLS Hedonic Regression ===")
print(ols.summary().tables[1])
print(f"\nAdjusted R²: {ols.rsquared_adj:.4f}")
# Marginal implicit prices at median house
median_price = df["price"].median()
print(f"\nMedian house price: ${median_price:,.0f}")
print(f"Implicit price of 1 extra sqft: ${ols.params['ln_sqft'] * median_price / df['sqft'].median():,.0f}")
print(f"CBD gradient: {ols.params['cbd_dist']*100:.1f}% per km")
# -----------------------------------------------------------------
# Bid-rent gradient visualization
# -----------------------------------------------------------------
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
# Price vs CBD distance
axes[0].scatter(df["cbd_dist"], df["ln_price"], alpha=0.3, s=5, c="steelblue")
d_line = np.linspace(0, 14, 100)
axes[0].plot(d_line, ols.params["const"] + ols.params["cbd_dist"] * d_line +
ols.params["ln_sqft"] * np.log(df["sqft"].median()) +
ols.params["bedrooms"] * 3 + ols.params["school_rating"] * 7,
"r-", lw=2, label="Bid-rent gradient")
axes[0].set_xlabel("Distance from CBD (km)")
axes[0].set_ylabel("ln(Price)")
axes[0].set_title("Bid-Rent Gradient")
axes[0].legend()
# Price heatmap on city grid
from scipy.interpolate import griddata
xi = np.linspace(0, 20, 100)
yi = np.linspace(0, 20, 100)
Xi, Yi = np.meshgrid(xi, yi)
Pi = griddata((df["x"], df["y"]), df["price"], (Xi, Yi), method="linear")
im = axes[1].contourf(Xi, Yi, Pi / 1e6, levels=20, cmap="RdYlGn_r")
axes[1].plot(10, 10, "k*", ms=15, label="CBD")
axes[1].set_title("Price Surface ($ million)")
plt.colorbar(im, ax=axes[1])
axes[1].legend()
# Residual spatial pattern
residuals = ols.resid
sc = axes[2].scatter(df["x"], df["y"], c=residuals, cmap="RdBu_r",
vmin=-0.3, vmax=0.3, s=10)
axes[2].plot(10, 10, "k*", ms=15)
axes[2].set_title("OLS Residuals (Spatial Pattern?)")
plt.colorbar(sc, ax=axes[2], label="Residual")
plt.tight_layout()
plt.savefig("hedonic_regression.png", dpi=150, bbox_inches="tight")
plt.close()
print("\nFigure saved: hedonic_regression.png")
Step 2: Repeat-Sales House Price Index
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
# -----------------------------------------------------------------
# Simulate repeat sales dataset
# 500 properties, each sold 2-3 times over 10 years (40 quarters)
# -----------------------------------------------------------------
np.random.seed(42)
n_properties = 500
n_quarters = 40
# True price index (simulates a housing boom-bust cycle)
true_index = np.zeros(n_quarters)
true_index[:20] = np.linspace(0, 0.8, 20) # 80% appreciation over 5 years
true_index[20:30] = np.linspace(0.8, 0.5, 10) # 30% crash
true_index[30:] = np.linspace(0.5, 0.7, 10) # partial recovery
records = []
for prop_id in range(n_properties):
# Two sales: first sale in 0-19, second in 10-39
t1 = np.random.randint(0, 25)
t2 = np.random.randint(t1 + 4, n_quarters)
base_ln_price = np.random.normal(12.5, 0.5) # property fixed effect
# Prices at each sale time
noise1 = np.random.normal(0, 0.08)
noise2 = np.random.normal(0, 0.08)
ln_p1 = base_ln_price + true_index[t1] + noise1
ln_p2 = base_ln_price + true_index[t2] + noise2
records.append({"prop_id": prop_id, "t": t1, "ln_price": ln_p1})
records.append({"prop_id": prop_id, "t": t2, "ln_price": ln_p2})
sales_df = pd.DataFrame(records).sort_values(["prop_id", "t"])
# -----------------------------------------------------------------
# Construct repeat-sales pairs
# -----------------------------------------------------------------
pairs = []
for prop_id, grp in sales_df.groupby("prop_id"):
grp = grp.reset_index(drop=True)
for i in range(len(grp) - 1):
pairs.append({
"prop_id": prop_id,
"t_sell": int(grp.loc[i+1, "t"]),
"t_buy": int(grp.loc[i, "t"]),
"ln_price_diff": grp.loc[i+1, "ln_price"] - grp.loc[i, "ln_price"]
})
pairs_df = pd.DataFrame(pairs)
print(f"Repeat-sale pairs: {len(pairs_df)}")
# -----------------------------------------------------------------
# Case-Shiller OLS repeat-sales regression
# Time dummies: omit period 0 as base
# -----------------------------------------------------------------
periods = list(range(1, n_quarters)) # periods 1..39
def build_repeat_sales_matrix(pairs_df, periods):
"""Build design matrix for repeat-sales regression.
Each row: ln(P_sell) - ln(P_buy) = sum of time dummies between buy and sell.
Returns X (design matrix) and y (price changes).
"""
n = len(pairs_df)
X = np.zeros((n, len(periods)))
y = pairs_df["ln_price_diff"].values.copy()
period_map = {p: i for i, p in enumerate(periods)}
for row_idx, row in enumerate(pairs_df.itertuples()):
for tau in range(row.t_buy + 1, row.t_sell + 1):
if tau in period_map:
X[row_idx, period_map[tau]] = 1
return X, y
X_rs, y_rs = build_repeat_sales_matrix(pairs_df, periods)
# OLS (no constant — natural normalization: index[0] = 0)
result_rs = sm.OLS(y_rs, X_rs).fit()
index_ols = np.concatenate([[0.0], result_rs.params]) # prepend 0 for period 0
print(f"R² of repeat-sales regression: {result_rs.rsquared:.4f}")
print(f"Peak index level: {index_ols.max():.3f} at period {np.argmax(index_ols)}")
# -----------------------------------------------------------------
# Visualization
# -----------------------------------------------------------------
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].plot(true_index, "b--", label="True Index", lw=2)
axes[0].plot(index_ols, "r-", label="OLS Repeat-Sales", lw=2)
axes[0].set_xlabel("Quarter")
axes[0].set_ylabel("Log Price Change (base=0)")
axes[0].set_title("Repeat-Sales House Price Index")
axes[0].legend()
axes[0].axhline(0, color="gray", lw=0.8)
# Distribution of holding period lengths
holding = pairs_df["t_sell"] - pairs_df["t_buy"]
axes[1].hist(holding, bins=20, color="steelblue", edgecolor="black")
axes[1].set_xlabel("Holding Period (quarters)")
axes[1].set_ylabel("Frequency")
axes[1].set_title("Distribution of Holding Periods")
plt.tight_layout()
plt.savefig("repeat_sales_index.png", dpi=150, bbox_inches="tight")
plt.close()
print("\nFigure saved: repeat_sales_index.png")
Step 3: Housing Affordability and Supply Analysis
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
# -----------------------------------------------------------------
# Simulate 20-city panel: 2000-2023 (annual)
# -----------------------------------------------------------------
np.random.seed(42)
n_cities = 20
n_years = 24
years = np.arange(2000, 2024)
city_names = [f"City_{i:02d}" for i in range(n_cities)]
# City characteristics
median_income_2000 = np.random.uniform(40000, 80000, n_cities)
# Income growth (varying by city)
income_growth = np.random.uniform(0.01, 0.04, n_cities)
# House price growth (some cities have higher appreciation)
price_growth = income_growth + np.random.uniform(-0.01, 0.04, n_cities)
# Housing supply elasticity (Saiz-style: constrained vs unconstrained)
supply_elast = np.random.uniform(0.5, 4.0, n_cities)
# Generate panel
records = []
for c in range(n_cities):
for t, yr in enumerate(years):
income = median_income_2000[c] * (1 + income_growth[c]) ** t
price = 200000 * (1 + price_growth[c]) ** t * (1 + np.random.normal(0, 0.02))
rent = price * 0.055 / 12 * (1 + np.random.normal(0, 0.02)) # ~5.5% cap rate
records.append({
"city": city_names[c],
"year": yr,
"median_income": income,
"median_price": price,
"median_rent": rent,
"supply_elasticity": supply_elast[c],
})
panel = pd.DataFrame(records)
# -----------------------------------------------------------------
# Affordability measures
# -----------------------------------------------------------------
INTEREST_RATE = 0.06 # 6% 30-year fixed
LOAN_TO_VALUE = 0.80
N_MONTHS = 360
def monthly_payment(price, rate=INTEREST_RATE, ltv=LOAN_TO_VALUE, n=N_MONTHS):
"""Compute monthly mortgage payment (principal + interest)."""
loan = price * ltv
r = rate / 12
return loan * r * (1 + r)**n / ((1 + r)**n - 1)
panel["monthly_pmt"] = panel["median_price"].apply(monthly_payment)
panel["annual_pmt"] = panel["monthly_pmt"] * 12
panel["PIR"] = panel["median_price"] / panel["median_income"]
panel["HAI"] = 100 * (panel["median_income"] / (panel["annual_pmt"] / 0.28))
panel["price_to_rent"] = panel["median_price"] / (panel["median_rent"] * 12)
panel["rent_burden"] = panel["median_rent"] / (panel["median_income"] / 12)
# Summary statistics
latest = panel[panel["year"] == 2023].copy()
print("=== 2023 Housing Affordability Summary ===")
print(latest[["city", "PIR", "HAI", "price_to_rent", "rent_burden"]].describe().round(2))
affordable = (latest["HAI"] >= 100).sum()
print(f"\nAffordable cities (HAI ≥ 100): {affordable}/{n_cities}")
# -----------------------------------------------------------------
# Supply elasticity and price growth relationship
# -----------------------------------------------------------------
price_growth_obs = latest.copy()
price_growth_obs["total_appreciation"] = (
panel.groupby("city").apply(
lambda g: (g.sort_values("year")["median_price"].iloc[-1] /
g.sort_values("year")["median_price"].iloc[0])
).values
)
corr, pval = stats.pearsonr(price_growth_obs["supply_elasticity"],
np.log(price_growth_obs["total_appreciation"]))
print(f"\nCorrelation(supply elasticity, ln price appreciation): r = {corr:.3f}, p = {pval:.3f}")
# -----------------------------------------------------------------
# Visualization
# -----------------------------------------------------------------
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# PIR over time by city
pivot_pir = panel.pivot(index="year", columns="city", values="PIR")
for col in pivot_pir.columns:
axes[0, 0].plot(pivot_pir.index, pivot_pir[col], alpha=0.4, lw=0.8, color="steelblue")
axes[0, 0].plot(pivot_pir.index, pivot_pir.mean(axis=1), "r-", lw=2, label="Mean PIR")
axes[0, 0].axhline(3.0, color="orange", ls="--", label="PIR=3 (affordable)")
axes[0, 0].set_title("Price-to-Income Ratio Over Time")
axes[0, 0].set_xlabel("Year"); axes[0, 0].set_ylabel("PIR")
axes[0, 0].legend()
# Housing Affordability Index
pivot_hai = panel.pivot(index="year", columns="city", values="HAI")
for col in pivot_hai.columns:
axes[0, 1].plot(pivot_hai.index, pivot_hai[col], alpha=0.4, lw=0.8, color="green")
axes[0, 1].plot(pivot_hai.index, pivot_hai.mean(axis=1), "k-", lw=2, label="Mean HAI")
axes[0, 1].axhline(100, color="red", ls="--", label="HAI=100 (threshold)")
axes[0, 1].set_title("Housing Affordability Index")
axes[0, 1].set_xlabel("Year"); axes[0, 1].set_ylabel("HAI")
axes[0, 1].legend()
# Supply elasticity vs price growth
axes[1, 0].scatter(price_growth_obs["supply_elasticity"],
price_growth_obs["total_appreciation"],
c="steelblue", edgecolors="k", s=60)
axes[1, 0].set_xlabel("Housing Supply Elasticity")
axes[1, 0].set_ylabel("Total Price Appreciation (2000-2023)")
axes[1, 0].set_title(f"Supply Constraint & Price Growth\n(r={corr:.2f}, p={pval:.3f})")
# Fit line
m, b = np.polyfit(price_growth_obs["supply_elasticity"],
price_growth_obs["total_appreciation"], 1)
xf = np.linspace(price_growth_obs["supply_elasticity"].min(),
price_growth_obs["supply_elasticity"].max(), 100)
axes[1, 0].plot(xf, m*xf + b, "r--")
# 2023 affordability distribution
axes[1, 1].hist(latest["PIR"], bins=10, color="coral", edgecolor="black", alpha=0.7,
label="PIR")
axes[1, 1].axvline(latest["PIR"].mean(), color="red", ls="--",
label=f"Mean PIR={latest['PIR'].mean():.1f}")
axes[1, 1].set_title("2023 PIR Distribution Across Cities")
axes[1, 1].set_xlabel("Price-to-Income Ratio")
axes[1, 1].legend()
plt.tight_layout()
plt.savefig("housing_affordability.png", dpi=150, bbox_inches="tight")
plt.close()
print("\nFigure saved: housing_affordability.png")
Advanced Usage
Spatial Lag Model for House Prices
import numpy as np
import pandas as pd
import statsmodels.api as sm
from scipy.spatial.distance import cdist
def spatial_weights_knn(coords, k=5):
"""Build row-standardized k-nearest neighbor spatial weight matrix.
Args:
coords: (n, 2) array of coordinates
k: number of neighbors
Returns:
W: (n, n) row-standardized weight matrix
"""
n = len(coords)
dist = cdist(coords, coords)
np.fill_diagonal(dist, np.inf)
W = np.zeros((n, n))
for i in range(n):
neighbors = np.argsort(dist[i])[:k]
W[i, neighbors] = 1
W = W / W.sum(axis=1, keepdims=True) # row standardize
return W
def spatial_lag_ols(y, X, W, max_iter=50, tol=1e-6):
"""Two-stage least squares estimator for spatial lag model.
Model: y = rho * W*y + X*beta + epsilon
Instruments: X, W*X for endogenous W*y term.
Args:
y: (n,) dependent variable
X: (n, k) exogenous regressors (with constant)
W: (n, n) spatial weight matrix
max_iter, tol: convergence criteria
Returns:
dict with rho, beta, residuals, fitted
"""
Wy = W @ y
WX = W @ X
# 2SLS: instrument Wy with WX
instruments = np.hstack([X, WX])
Z_hat = instruments @ np.linalg.pinv(instruments) @ np.column_stack([X, Wy])
beta_2sls = np.linalg.lstsq(Z_hat, y, rcond=None)[0]
rho = beta_2sls[-1]
beta = beta_2sls[:-1]
fitted = rho * Wy + X @ beta
residuals = y - fitted
return {"rho": rho, "beta": beta, "residuals": residuals, "fitted": fitted}
np.random.seed(42)
n_h = 400
coords_h = np.random.uniform(0, 20, (n_h, 2))
W_h = spatial_weights_knn(coords_h, k=8)
cbd_d = np.sqrt((coords_h[:, 0] - 10)**2 + (coords_h[:, 1] - 10)**2)
sqft_h = np.random.lognormal(7, 0.3, n_h)
# True DGP includes spatial lag
rho_true = 0.4
X_h = np.column_stack([np.ones(n_h), np.log(sqft_h), cbd_d])
eps_h = np.random.normal(0, 0.1, n_h)
ln_p_h = np.linalg.solve(np.eye(n_h) - rho_true * W_h,
X_h @ np.array([11, 0.7, -0.04]) + eps_h)
result_sl = spatial_lag_ols(ln_p_h, X_h, W_h)
print("=== Spatial Lag Model (2SLS) ===")
print(f"Spatial autoregressive parameter ρ = {result_sl['rho']:.4f} (true: {rho_true})")
print(f"Coefficients: {result_sl['beta'].round(4)}")
Capitalization Rate and Yield Analysis
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
def compute_yields(purchase_price, annual_rent, operating_expenses_ratio=0.40,
appreciation_rate=0.03, hold_years=5, discount_rate=0.08):
"""Compute key real estate investment metrics.
Args:
purchase_price: acquisition price
annual_rent: gross annual rent
operating_expenses_ratio: fraction of rent for expenses (NOI/rent = 1-OER)
appreciation_rate: annual property appreciation
hold_years: investment horizon in years
discount_rate: required return for NPV
Returns:
dict of metrics
"""
noi = annual_rent * (1 - operating_expenses_ratio)
cap_rate = noi / purchase_price
gross_yield = annual_rent / purchase_price
net_yield = noi / purchase_price
# Cash flows assuming NOI grows at appreciation rate
cash_flows = [-purchase_price]
for t in range(1, hold_years + 1):
noi_t = noi * (1 + appreciation_rate) ** (t - 1)
if t < hold_years:
cash_flows.append(noi_t)
else:
# Terminal sale: resale price based on terminal cap rate
terminal_noi = noi * (1 + appreciation_rate) ** t
resale_price = terminal_noi / cap_rate
cash_flows.append(noi_t + resale_price)
# NPV and IRR
npv = sum(cf / (1 + discount_rate)**t for t, cf in enumerate(cash_flows))
# IRR via Newton's method
irr = discount_rate
for _ in range(200):
npv_irr = sum(cf / (1 + irr)**t for t, cf in enumerate(cash_flows))
dnpv = sum(-t * cf / (1 + irr)**(t+1) for t, cf in enumerate(cash_flows))
if abs(dnpv) < 1e-10:
break
irr -= npv_irr / dnpv
return {
"cap_rate": cap_rate,
"gross_yield": gross_yield,
"net_yield": net_yield,
"npv": npv,
"irr": irr,
"price_to_rent": purchase_price / annual_rent,
}
metrics = compute_yields(500000, 24000, operating_expenses_ratio=0.35, discount_rate=0.07)
print("=== Investment Property Metrics ===")
for k, v in metrics.items():
if k in ("cap_rate", "gross_yield", "net_yield", "irr"):
print(f" {k}: {v*100:.2f}%")
elif k == "npv":
print(f" {k}: ${v:,.0f}")
else:
print(f" {k}: {v:.1f}")
Troubleshooting
| Problem | Cause | Fix |
|---|---|---|
| Negative fitted prices in linear model | Outliers pulling fit | Use log-linear specification |
| Multicollinearity (sqft ↔ bedrooms) | Structural correlation | Check VIF; consider sqft per room instead |
| Spatial autocorrelation in residuals | Omitted spatial effects | Add spatial lag or spatial error term |
| Repeat-sales singular matrix | Too few observations per period | Aggregate to quarterly from monthly |
| IRR computation diverges | No sign change in cash flows | Check if investment is ever profitable |
| HAI > 200 (implausibly high) | Low interest rate scenario | Verify mortgage rate inputs |
External Resources
- Rosen, S. (1974). Hedonic prices and implicit markets. Journal of Political Economy, 82(1).
- Case, K. E., & Shiller, R. J. (1989). The efficiency of the market for single-family homes. AER, 79(1).
- Saiz, A. (2010). The geographic determinants of housing supply. QJE, 125(3).
- Anselin, L. (1988). Spatial Econometrics. Kluwer Academic.
- S&P Case-Shiller Index methodology
- Zillow Research Data — ZHVI, ZRI
Examples
Example 1: Neighborhood Price Gradient by Amenity
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
np.random.seed(42)
n = 800
x = np.random.uniform(0, 15, n)
y = np.random.uniform(0, 15, n)
# Two amenities: park at (3, 8) and school at (12, 5)
park_dist = np.sqrt((x - 3)**2 + (y - 8)**2)
school_dist = np.sqrt((x - 12)**2 + (y - 5)**2)
cbd_dist = np.sqrt((x - 7.5)**2 + (y - 7.5)**2)
sqft = np.random.lognormal(6.8, 0.35, n)
ln_p = (11.5
+ 0.65 * np.log(sqft)
- 0.035 * cbd_dist
- 0.045 * park_dist # park premium: $45k per km closer
- 0.020 * school_dist
+ np.random.normal(0, 0.10, n))
df = pd.DataFrame({"ln_p": ln_p, "ln_sqft": np.log(sqft),
"cbd_dist": cbd_dist, "park_dist": park_dist,
"school_dist": school_dist})
X = sm.add_constant(df[["ln_sqft", "cbd_dist", "park_dist", "school_dist"]])
res = sm.OLS(df["ln_p"], X).fit(cov_type="HC3")
med_p = np.exp(df["ln_p"].mean())
print("=== Amenity Valuation ===")
print(f"Park premium per 1km closer: ${abs(res.params['park_dist']) * med_p:,.0f}")
print(f"School proximity: ${abs(res.params['school_dist']) * med_p:,.0f} per km")
print(f"CBD gradient: {res.params['cbd_dist']*100:.2f}% per km")
Example 2: Housing Market Bubble Detection (Price-to-Rent Ratio)
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
np.random.seed(2024)
quarters = pd.date_range("2000Q1", "2023Q4", freq="QE")
n_q = len(quarters)
# Simulate price and rent indices (rent grows more smoothly)
rent_idx = 100 * (1.025 ** (np.arange(n_q) / 4)) # 2.5%/yr
# Price: boom-bust cycle
price_growth = np.concatenate([
np.linspace(0, 0.10, 24), # Q1 2000 - Q4 2005 boom
np.linspace(0.10, 0.02, 12), # Q1 2006 - Q4 2008 peak/bust
np.linspace(0.02, -0.05, 8), # Q1 2009 - Q4 2010 trough
np.linspace(-0.05, 0.05, 24), # recovery
np.linspace(0.05, 0.12, n_q - 68), # 2017-2023 new boom
])[:n_q]
price_idx = 100 * np.cumprod(1 + price_growth / 4)
ptr = price_idx / rent_idx * (rent_idx[0] / price_idx[0]) # normalized P/R ratio
# Hodrick-Prescott trend (approximate with polynomial)
from numpy.polynomial import polynomial as P
t = np.arange(n_q)
coeffs = P.polyfit(t, ptr, 3)
ptr_trend = P.polyval(t, coeffs)
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(quarters, ptr, "steelblue", label="Price-to-Rent Ratio")
ax.plot(quarters, ptr_trend, "r--", label="Trend (cubic)")
ax.axhline(1.0, color="gray", ls=":", label="Baseline")
overvalued = ptr > ptr_trend * 1.10
ax.fill_between(quarters, ptr, ptr_trend,
where=overvalued, alpha=0.3, color="red", label="Overvalued >10%")
ax.set_title("Housing Market Bubble Detection via Price-to-Rent Ratio")
ax.set_ylabel("Normalized P/R Ratio")
ax.legend()
plt.tight_layout()
plt.savefig("bubble_detection.png", dpi=150, bbox_inches="tight")
plt.close()
print(f"Overvalued quarters: {overvalued.sum()}/{n_q}")