TeX 55.1%
Python 44.2%
1# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca2"""Synthetic Quebec-like housing markets with known quality-constant latent price paths (protocol §12).34A market has R regions × M MRCs × municipalities of heterogeneous size (expected monthly sales from a log-normal5mimicking results/profile/municipality_size_bands.csv). True log index of municipality u at month t:6 L_{u,t} = μ_t + δ^R_{r(u),t} + δ^M_{m(u),t} + δ^U_{u,t}7with μ a random walk with drift and optional structural break / turning point, and nested random-walk deviations8whose innovation sd differs across municipalities. Prices:9 log P_i = L_{u(i),t(i)} + f_t(X_i) + α_cell(i) + ε_i, ε ~ Student-t(4) scaled (+ contamination)10f is nonlinear in floor area with an optional drifting elasticity; α_cell is spatially correlated.11Composition shifts: the distribution of sold X and of sold cells can trend over time (quality-up / quality-down,12location drift). Sales counts are Poisson with June–July seasonality. Floor area missing MCAR or time-increasing.13Repeat sales: a fraction of parcels resell. Registration lags: exponential days (for real-time experiments).14The output has the same columns as the real analytic layer so every estimator runs unchanged.15"""16from __future__ import annotations1718from dataclasses import dataclass, field1920import numpy as np21import polars as pl222324@dataclass25class Scenario:26 name: str = "baseline"27 T: int = 6028 n_regions: int = 429 mrc_per_region: int = 330 mun_per_mrc: int = 531 size_logmean: float = 1.2 # log expected monthly sales per municipality (median ≈ 3.3)32 size_logsd: float = 1.333 drift: float = 0.004 # provincial monthly drift34 tau_prov: float = 0.01235 tau_region: float = 0.00636 tau_mrc: float = 0.00437 tau_mun: float = 0.00638 tau_mun_hetero: bool = True # half the municipalities have tau_mun/439 break_month: int | None = None # turning point: drift sign flips40 break_drift: float = -0.00841 sigma_eps: float = 0.1642 t_df: float = 4.043 contamination: float = 0.0144 beta_fa: float = 0.5545 beta_fa_drift: float = 0.0 # change in floor-area elasticity over the sample46 quality_trend: float = 0.0 # monthly drift of mean log floor area of SOLD homes (composition)47 location_trend: float = 0.0 # monthly drift of mean cell quality of SOLD homes48 fa_missing: float = 0.0349 fa_missing_trend: float = 0.0 # additional missingness at the end of the sample50 resale_prob: float = 0.1251 seasonal_amp: float = 0.652 reg_lag_mean_days: float = 10.053 seed: int = 054 scale_sales: float = 1.0 # multiply all expected sales (sample-size experiments)555657@dataclass58class SimOutput:59 df: pl.DataFrame60 truth_mun: np.ndarray # (n_mun, T) true log index61 truth_region: np.ndarray # (n_regions, T) stock-weighted62 truth_mrc: np.ndarray63 truth_prov: np.ndarray # (T,)64 mun_codes: list[str]65 mun_region: np.ndarray66 mun_mrc: np.ndarray67 stock: np.ndarray # housing stock weights per municipality68 months: list[str]69 scenario: Scenario707172def _months(T: int, start: str = "2021-01") -> list[str]:73 y, m = int(start[:4]), int(start[5:7]); out = []74 for _ in range(T):75 out.append(f"{y:04d}-{m:02d}"); m += 176 if m == 13: y, m = y + 1, 177 return out787980def simulate(sc: Scenario) -> SimOutput:81 rng = np.random.default_rng(sc.seed)82 T = sc.T; months = _months(T)83 n_mrc = sc.n_regions * sc.mrc_per_region; n_mun = n_mrc * sc.mun_per_mrc84 mun_mrc = np.repeat(np.arange(n_mrc), sc.mun_per_mrc); mun_region = mun_mrc // sc.mrc_per_region85 lam = np.exp(rng.normal(sc.size_logmean, sc.size_logsd, n_mun)) * sc.scale_sales86 stock = lam * rng.uniform(80, 160, n_mun) # stock ∝ turnover with noise87 # latent paths88 drift = np.full(T, sc.drift)89 if sc.break_month is not None:90 drift[sc.break_month:] = sc.break_drift91 mu = np.cumsum(drift + rng.normal(0, sc.tau_prov, T)); mu -= mu[0]92 dR = np.cumsum(rng.normal(0, sc.tau_region, (sc.n_regions, T)), axis=1); dR -= dR[:, [0]]93 dM = np.cumsum(rng.normal(0, sc.tau_mrc, (n_mrc, T)), axis=1); dM -= dM[:, [0]]94 tau_u = np.full(n_mun, sc.tau_mun)95 if sc.tau_mun_hetero:96 tau_u[rng.random(n_mun) < 0.5] = sc.tau_mun / 497 dU = np.cumsum(rng.normal(0, 1, (n_mun, T)) * tau_u[:, None], axis=1); dU -= dU[:, [0]]98 L = mu[None, :] + dR[mun_region] + dM[mun_mrc] + dU99 # geography: municipality centres clustered by region, cells inside100 reg_c = rng.uniform(0, 400_000, (sc.n_regions, 2))101 mrc_c = reg_c[np.arange(n_mrc) // sc.mrc_per_region] + rng.normal(0, 40_000, (n_mrc, 2))102 mun_c = mrc_c[mun_mrc] + rng.normal(0, 12_000, (n_mun, 2))103 mun_level = rng.normal(0, 0.25, n_mun) # price level of the municipality (absorbed by location FE)104 # sales105 season = 1 + sc.seasonal_amp * np.cos(2 * np.pi * (np.arange(T) % 12 - 5) / 12)106 rows = []107 parcel_counter = 0108 for u in range(n_mun):109 for t in range(T):110 n = rng.poisson(lam[u] * season[t])111 if n == 0:112 continue113 # composition drift of sold homes114 lfa = rng.normal(np.log(120) + 0.15 * mun_level[u] + sc.quality_trend * t, 0.35, n)115 age = np.clip(rng.gamma(2.0, 20, n), 0, 150)116 lla = rng.normal(np.log(600), 0.6, n)117 storeys = rng.choice([1, 2], n, p=[0.6, 0.4])118 # within-municipality location: cell offset with spatially correlated quality; location drift = mean cell quality trends119 theta = rng.uniform(0, 2 * np.pi, n); rad = rng.gamma(2, 1200, n)120 x = mun_c[u, 0] + rad * np.cos(theta); y = mun_c[u, 1] + rad * np.sin(theta)121 cell_q = 0.08 * np.sin(x / 1500) * np.cos(y / 1500) + 0.06 * (rad < 1500) + sc.location_trend * t * rng.uniform(0, 1, n)122 beta_fa_t = sc.beta_fa + sc.beta_fa_drift * t / max(T - 1, 1)123 fX = beta_fa_t * (lfa - np.log(120)) + 0.08 * (lfa - np.log(120)) ** 2 - 0.004 * age + 0.00002 * age ** 2 + 0.09 * (lla - np.log(600)) + 0.05 * (storeys - 1)124 eps = rng.standard_t(sc.t_df, n) * sc.sigma_eps / np.sqrt(sc.t_df / (sc.t_df - 2))125 cont = rng.random(n) < sc.contamination126 eps[cont] += rng.choice([-1.0, 1.0], cont.sum()) * rng.uniform(0.6, 1.2, cont.sum())127 lp = 12.6 + mun_level[u] + L[u, t] + fX + cell_q + eps128 pid = parcel_counter + np.arange(n); parcel_counter += n129 rows.append(pl.DataFrame({130 "mun_i": u, "t": t, "log_price": lp, "lfa_true": lfa, "age": age, "lla_true": lla, "storeys": storeys, "x_m": x, "y_m": y, "parcel_i": pid,131 "year_built_offset": age,132 }))133 df = pl.concat(rows)134 # repeat sales: a share of sales are resales of an earlier parcel (same X, same cell, new time)135 n_rows = df.height136 is_resale = rng.random(n_rows) < sc.resale_prob137 df = df.with_row_index("_i")138 # pick an earlier parcel in the same municipality for resales139 first = df.filter(~pl.Series(is_resale)).select(["_i", "mun_i", "t", "parcel_i", "lfa_true", "age", "lla_true", "storeys", "x_m", "y_m"])140 res = df.filter(pl.Series(is_resale)).select(["_i", "mun_i", "t", "log_price"])141 cand = first.join(res, on="mun_i", how="inner", suffix="_r").filter(pl.col("t_r") - pl.col("t") >= 6)142 if cand.height: # one random earlier parcel per resale (not always the earliest)143 cand = cand.with_columns(pl.Series("_rk", rng.random(cand.height))).sort("_rk").group_by("_i_r").agg(pl.all().first()).drop("_rk")144 if cand.height:145 # replace resale rows' parcel/X by the matched earlier sale (price keeps the new period's latent level)146 rep = cand.select(pl.col("_i_r").alias("_i"), pl.col("parcel_i").alias("parcel_new"), pl.col("lfa_true").alias("lfa_new"), pl.col("age").alias("age_new"),147 pl.col("lla_true").alias("lla_new"), pl.col("storeys").alias("st_new"), pl.col("x_m").alias("x_new"), pl.col("y_m").alias("y_new"))148 df = df.join(rep, on="_i", how="left")149 # price adjustment so the resale keeps its own latent level: recompute hedonic part difference approximately via lfa change150 beta = sc.beta_fa151 df = df.with_columns(152 pl.when(pl.col("parcel_new").is_not_null()).then(pl.col("log_price") + beta * (pl.col("lfa_new") - pl.col("lfa_true"))).otherwise(pl.col("log_price")).alias("log_price"),153 pl.coalesce(pl.col("parcel_new"), pl.col("parcel_i")).alias("parcel_i"),154 pl.coalesce(pl.col("lfa_new"), pl.col("lfa_true")).alias("lfa_true"), pl.coalesce(pl.col("age_new"), pl.col("age")).alias("age"),155 pl.coalesce(pl.col("lla_new"), pl.col("lla_true")).alias("lla_true"), pl.coalesce(pl.col("st_new"), pl.col("storeys")).alias("storeys"),156 pl.coalesce(pl.col("x_new"), pl.col("x_m")).alias("x_m"), pl.coalesce(pl.col("y_new"), pl.col("y_m")).alias("y_m"),157 ).drop(["parcel_new", "lfa_new", "age_new", "lla_new", "st_new", "x_new", "y_new"])158 # missing floor area (MCAR + trend)159 p_miss = sc.fa_missing + sc.fa_missing_trend * (df["t"].to_numpy() / max(T - 1, 1))160 miss = rng.random(df.height) < p_miss161 mun_codes = [f"M{u:03d}" for u in range(n_mun)]162 mrc_codes = [f"R{m:02d}" for m in range(n_mrc)]163 reg_codes = [f"{r:02d}" for r in range(sc.n_regions)]164 lag_days = rng.exponential(sc.reg_lag_mean_days, df.height)165 df = df.with_columns(166 pl.Series("month", [months[t] for t in df["t"].to_list()]),167 pl.Series("mun_code", [mun_codes[u] for u in df["mun_i"].to_list()]),168 pl.Series("mrc_code", [mrc_codes[mun_mrc[u]] for u in df["mun_i"].to_list()]),169 pl.Series("region_code", [reg_codes[mun_region[u]] for u in df["mun_i"].to_list()]),170 pl.Series("floor_area_m2", np.where(miss, np.nan, np.exp(df["lfa_true"].to_numpy()))),171 pl.col("lla_true").exp().alias("land_area_m2"),172 pl.Series("reg_lag_days", lag_days),173 ).with_columns(174 pl.col("month").str.to_date("%Y-%m-%d", strict=False).alias("_d"),175 )176 # dates: uniform day within month177 day = rng.integers(1, 28, df.height)178 df = df.with_columns(pl.Series("_day", day)).with_columns(179 (pl.col("month") + "-" + pl.col("_day").cast(pl.Utf8).str.zfill(2)).str.to_date("%Y-%m-%d").alias("date"),180 ).with_columns(181 pl.col("date").dt.year().alias("year"), pl.col("log_price").exp().alias("price"),182 pl.lit("house").alias("property_class"), pl.lit("detached").alias("subtype"), pl.lit("De plain-pied").alias("genre"),183 pl.lit(1).alias("n_units"), pl.col("storeys").cast(pl.Float64).alias("storeys"),184 pl.col("parcel_i").cast(pl.Utf8).alias("parcel_id"), pl.col("_i").cast(pl.Utf8).alias("tx_id"),185 pl.when(pl.Series(miss)).then(None).otherwise(pl.col("floor_area_m2")).alias("floor_area_m2"),186 (pl.col("date") + pl.duration(days=pl.col("reg_lag_days").cast(pl.Int64))).alias("available_date"),187 )188 for lvl, size in (("grid500", 500), ("grid1k", 1000), ("grid2k", 2000), ("grid5k", 5000)):189 df = df.with_columns((pl.col("x_m") // size).cast(pl.Int64).cast(pl.Utf8).add("_").add((pl.col("y_m") // size).cast(pl.Int64).cast(pl.Utf8)).alias(lvl))190 df = df.drop(["_i", "_d", "_day", "mun_i", "t", "parcel_i", "lfa_true", "lla_true", "year_built_offset"])191 # truths: stock-weighted aggregates192 w = stock / stock.sum()193 truth_prov = (w[:, None] * L).sum(axis=0)194 truth_region = np.vstack([(w[mun_region == r][:, None] * L[mun_region == r]).sum(axis=0) / w[mun_region == r].sum() for r in range(sc.n_regions)])195 truth_mrc = np.vstack([(w[mun_mrc == m][:, None] * L[mun_mrc == m]).sum(axis=0) / w[mun_mrc == m].sum() for m in range(n_mrc)])196 return SimOutput(df, L, truth_region, truth_mrc, truth_prov, mun_codes, mun_region, mun_mrc, stock, months, sc)197198199SCENARIOS: dict[str, Scenario] = {200 "baseline": Scenario(name="baseline"),201 "turning_point": Scenario(name="turning_point", break_month=36, drift=0.006, break_drift=-0.008),202 "flat_quality_up": Scenario(name="flat_quality_up", drift=0.0, tau_prov=0.006, quality_trend=0.004, location_trend=0.0008),203 "rising_quality_down": Scenario(name="rising_quality_down", drift=0.008, quality_trend=-0.004, location_trend=-0.0008),204 "drifting_hedonics": Scenario(name="drifting_hedonics", beta_fa_drift=0.15, fa_missing_trend=0.08),205 "heavy_tails": Scenario(name="heavy_tails", t_df=2.5, contamination=0.04, sigma_eps=0.2),206 "sparse": Scenario(name="sparse", scale_sales=0.35),207 "no_local_dynamics": Scenario(name="no_local_dynamics", tau_mun=0.0005, tau_mrc=0.0005, tau_region=0.001), # counter-example: pooling should win big208 "idiosyncratic_towns": Scenario(name="idiosyncratic_towns", tau_mun=0.02, tau_mun_hetero=False, tau_region=0.002), # counter-example: pooling hurts209}210