Real Estate Market Analysis

import pandas as pd
import numpy as np
from scipy import stats

df = pd.DataFrame({
    "prop_id":      [f"P{i:03d}" for i in range(1,21)],
    "price":        [485000,218000,742000,195000,368000,
                     312000,895000,155000,428000,284000,
                     562000,176000,318000,680000,245000,
                     415000,720000,198000,342000,524000],
    "sqft":         [1240,680,1820,610,980,
                     840,2100,520,1100,760,
                     1380,590,870,1680,650,
                     1050,1940,620,920,1280],
    "bedrooms":     [3,2,4,1,3,2,5,1,3,2,4,2,3,4,2,3,4,2,3,3],
    "bathrooms":    [2,1,3,1,2,1,4,1,2,1,3,1,2,3,1,2,3,1,2,2],
    "age_years":    [12,34,8,48,22,18,5,62,14,28,10,44,20,7,38,16,6,52,24,11],
    "prop_type":    ["Semi-detached","Flat","Detached","Flat","Terraced",
                     "Semi-detached","Detached","Flat","Semi-detached","Flat",
                     "Detached","Flat","Terraced","Detached","Flat",
                     "Semi-detached","Detached","Flat","Terraced","Semi-detached"],
    "neighbourhood":["Kensington","Hackney","Kensington","Hackney","Islington",
                      "Islington","Kensington","Hackney","Kensington","Hackney",
                      "Kensington","Hackney","Islington","Kensington","Hackney",
                      "Islington","Kensington","Hackney","Islington","Kensington"],
    "condition":    [4,3,5,2,4,3,5,2,4,3,5,2,4,5,3,4,5,2,3,4]
})

# ── Log-transform price (dependent variable) ─────────────────────────────────
df["log_price"]   = np.log(df["price"])
df["price_per_sf"]= (df["price"] / df["sqft"]).round(1)

# ── Dummy encoding: drop_first=True avoids perfect multicollinearity ──────────
dummies = pd.get_dummies(df[["prop_type","neighbourhood"]], drop_first=True)
dummies = dummies.astype(int)   # convert bool to int for matrix ops

print("Dataset shape:", df.shape)
print(f"\nPrice summary:")
print(df["price"].describe().apply(lambda x: f"£{x:,.0f}").to_string())
print(f"\nDummy columns created: {list(dummies.columns)}")
print(f"\nPrice per sq ft by neighbourhood:")
print(df.groupby("neighbourhood")["price_per_sf"].agg(["mean","min","max"]).round(0).to_string())
print(f"\nCorrelations with log_price:")
for col in ["sqft","bedrooms","bathrooms","age_years","condition"]:
    r = df[col].corr(df["log_price"]).round(3)
    print(f"  {col:<14} r = {r:+.3f}")

# ── Build feature matrix X ────────────────────────────────────────────────────
continuous = df[["sqft","bedrooms","bathrooms","age_years","condition"]].copy()
continuous["sqft_100"] = continuous["sqft"] / 100   # scale: coefficient = % per 100 sqft

X_raw = pd.concat([
    continuous[["sqft_100","bedrooms","bathrooms","age_years","condition"]],
    dummies
], axis=1)

# Add intercept column (column of ones)
X = np.column_stack([np.ones(len(X_raw)), X_raw.values])
y = df["log_price"].values
feature_names = ["intercept"] + list(X_raw.columns)

# ── OLS solution via numpy least squares ──────────────────────────────────────
# beta = (X'X)^{-1} X'y  — numpy lstsq solves this efficiently
beta, residuals_ss, rank, sv = np.linalg.lstsq(X, y, rcond=None)

# ── Residuals and model fit ───────────────────────────────────────────────────
y_hat      = X @ beta
resid      = y - y_hat
n, k       = len(y), len(beta)
dof        = n - k              # degrees of freedom
sigma2     = resid @ resid / dof   # residual variance

# R-squared
ss_res = resid @ resid
ss_tot = ((y - y.mean()) ** 2).sum()
r2     = 1 - ss_res / ss_tot
r2_adj = 1 - (1-r2) * (n-1) / (n-k-1)

# ── Standard errors, t-statistics, p-values ───────────────────────────────────
XtX_inv   = np.linalg.inv(X.T @ X)
se        = np.sqrt(sigma2 * np.diag(XtX_inv))
t_stat    = beta / se
p_val     = 2 * (1 - stats.t.cdf(np.abs(t_stat), df=dof))
ci_lo     = beta - 1.96 * se
ci_hi     = beta + 1.96 * se

print(f"Hedonic Pricing Model — OLS Results")
print(f"R² = {r2:.4f}  |  Adjusted R² = {r2_adj:.4f}  |  n = {n}  |  df = {dof}")
print(f"\n{'Feature':<28} {'Coeff':>8} {'SE':>7} {'t-stat':>8} {'p-val':>8} {'95% CI':>22} {'Sig'}")
print("─" * 90)
for i, feat in enumerate(feature_names):
    sig = "***" if p_val[i]<0.001 else "**" if p_val[i]<0.01 else "*" if p_val[i]<0.05 else ""
    pct = f"({(np.exp(beta[i])-1)*100:+.1f}%)" if feat != "intercept" else ""
    print(f"  {feat:<26} {beta[i]:>8.4f} {se[i]:>7.4f} {t_stat[i]:>8.2f} {p_val[i]:>8.4f}  [{ci_lo[i]:.3f}, {ci_hi[i]:.3f}]  {sig} {pct}")

# ── Predicted prices and residuals ────────────────────────────────────────────
df["log_price_pred"] = y_hat
df["price_pred"]     = np.exp(y_hat).round(0)
df["log_resid"]      = resid
df["pct_resid"]      = ((df["price"] - df["price_pred"]) / df["price_pred"] * 100).round(1)
df["upside_gbp"]     = (df["price_pred"] - df["price"]).round(0)   # positive = undervalued

# ── Full residual table ────────────────────────────────────────────────────────
result = df[["prop_id","neighbourhood","prop_type","sqft","price",
             "price_pred","pct_resid","upside_gbp"]].copy()
result = result.sort_values("pct_resid")  # most undervalued first

print("Property valuation — sorted by mispricing (most undervalued first):")
print(result.to_string(index=False))

# ── Top 3 acquisition targets ─────────────────────────────────────────────────
targets = result[result["pct_resid"] < 0].head(3)
print(f"\nTop 3 acquisition targets (largest upside):")
for _, row in targets.iterrows():
    print(f"  {row['prop_id']}  {row['neighbourhood']:<12} {row['prop_type']:<16}"
          f"  Actual: £{row['price']:,.0f}  Predicted: £{row['price_pred']:,.0f}"
          f"  Upside: £{row['upside_gbp']:,.0f} ({-row['pct_resid']:.1f}% below model)")

# ── Residual normality test ────────────────────────────────────────────────────
_, p_shapiro = stats.shapiro(resid)
print(f"\nResidual diagnostics:")
print(f"  Mean residual:      {resid.mean():+.6f} (should be ~0)")
print(f"  Std residual:       {resid.std():.4f}")
print(f"  Shapiro-Wilk p:     {p_shapiro:.4f}  ({'Normal ✓' if p_shapiro>0.05 else 'Non-normal'})")

# ── Standardised coefficients (exclude intercept and dummies) ─────────────────
cont_features = ["sqft_100","bedrooms","bathrooms","age_years","condition"]
cont_idx      = [feature_names.index(f) for f in cont_features]

X_std_df = X_raw[cont_features].copy()
X_stds   = X_std_df.std().values
beta_std = beta[cont_idx] * X_stds   # standardised: β × std(X)

print("Standardised feature importance (continuous variables):")
imp_df = pd.DataFrame({
    "feature":   cont_features,
    "coeff":     beta[cont_idx].round(4),
    "std_coeff": np.abs(beta_std).round(4),
    "pct_premium": ((np.exp(beta[cont_idx])-1)*100).round(1)
}).sort_values("std_coeff", ascending=False)
print(imp_df.to_string(index=False))

# ── Value-add scenarios ────────────────────────────────────────────────────────
# Renovation uplift: condition improvement from 2 to 5 (+3 condition points)
condition_coeff = beta[feature_names.index("condition")]
renovation_uplift_pct = (np.exp(condition_coeff * 3) - 1) * 100

# On a Hackney flat (mean price £205k) vs Kensington detached (mean price £720k)
mean_hackney    = df[df["neighbourhood"]=="Hackney"]["price"].mean()
mean_kensington = df[(df["neighbourhood"]=="Kensington") &
                     (df["prop_type"]=="Detached")]["price"].mean()

print(f"\nValue-add analysis — full renovation (condition 2→5, +3 points):")
print(f"  Renovation uplift:     {renovation_uplift_pct:.1f}%")
print(f"  On Hackney flat (~£{mean_hackney:,.0f}):   +£{mean_hackney*renovation_uplift_pct/100:,.0f}")
print(f"  On Kensington detached (~£{mean_kensington:,.0f}): +£{mean_kensington*renovation_uplift_pct/100:,.0f}")

# Age discount: 30-year-old property vs new-build
age_coeff = beta[feature_names.index("age_years")]
age_discount_pct = (np.exp(age_coeff * 30) - 1) * 100
print(f"\nAge discount — 30 years older:")
print(f"  Implied discount:      {age_discount_pct:.1f}%  (vs new-build)")

# ── Leverage and Cook's distance (influential observations) ───────────────────
H        = X @ np.linalg.inv(X.T @ X) @ X.T   # hat matrix
leverage = np.diag(H)                          # h_ii: influence of each obs
std_resid= resid / (sigma2**0.5 * np.sqrt(1 - leverage))  # studentised residuals

# Cook's distance: combined leverage × residual measure
k_params  = len(beta)
cooks_d   = (std_resid**2 * leverage) / (k_params * (1 - leverage))
df["leverage"]  = leverage.round(4)
df["std_resid"] = std_resid.round(4)
df["cooks_d"]   = cooks_d.round(4)

print("Influential observation diagnostics:")
print(df[["prop_id","price","pct_resid","leverage","std_resid","cooks_d"]]
      .sort_values("cooks_d", ascending=False).head(5).to_string(index=False))

# ── VIF: variance inflation factor for multicollinearity ─────────────────────
print("\nVariance Inflation Factors (continuous features):")
for i, feat in enumerate(cont_features):
    # VIF = 1 / (1 - R² of regressing this feature on all others)
    X_others = np.delete(X_raw[cont_features].values, i, axis=1)
    X_o = np.column_stack([np.ones(len(X_others)), X_others])
    b_o, _, _, _ = np.linalg.lstsq(X_o, X_raw[feat].values, rcond=None)
    y_hat_o = X_o @ b_o
    ss_r = ((X_raw[feat].values - y_hat_o)**2).sum()
    ss_t = ((X_raw[feat].values - X_raw[feat].values.mean())**2).sum()
    r2_o = 1 - ss_r/ss_t
    vif  = 1 / (1 - r2_o) if r2_o < 1 else np.inf
    flag = " ← high (>5)" if vif > 5 else ""
    print(f"  {feat:<14}  VIF = {vif:.2f}{flag}")

# ── Acquisition brief ─────────────────────────────────────────────────────────
print("\n" + "="*60)
print("INVESTMENT COMMITTEE BRIEF — TOP 3 ACQUISITION TARGETS")
print("="*60)
targets_brief = df.nsmallest(3, "pct_resid")[
    ["prop_id","neighbourhood","prop_type","sqft","condition","price",
     "price_pred","pct_resid","upside_gbp"]
]
for _, r in targets_brief.iterrows():
    print(f"\n  {r['prop_id']} — {r['neighbourhood']} {r['prop_type']}, {r['sqft']} sqft, Condition {r['condition']}")
    print(f"    Asking price:   £{r['price']:>10,.0f}")
    print(f"    Model value:    £{r['price_pred']:>10,.0f}")
    print(f"    Discount:            {-r['pct_resid']:.1f}% below model")
    print(f"    Mispricing upside:  £{r['upside_gbp']:>10,.0f}")
    reno = r['price_pred'] * (np.exp(condition_coeff*(5-r['condition']))-1)
    print(f"    Renovation upside:  £{reno:>10,.0f}  (condition {r['condition']:.0f}→5)")