TeX 62.3%
Python 37.6%
1# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca2"""Forecasting datasets built strictly from information available at each origin.34Task A (property level): for origin o, the training set is every headline sale dated ≤ o with features computed from5sales ≤ o; the test set is the sales of months o+1 … o+12, each labelled by its horizon h = months after o. Market6features for a test sale are those of its municipality *at the origin* (not at the sale date).7Task B (cell level): for origin o and cell g, predictors are the real-time SHARP-HPI vintage published for o (WP228index_vintages, smoothed mode), volumes and dispersion up to o; target = final smoothed log index change over h.9"""10from __future__ import annotations1112import sys1314import numpy as np15import polars as pl16from scipy.spatial import cKDTree1718from wp23 import config as C1920sys.path.insert(0, str(C.WP22_SRC))21from sharp_hpi.hedonic import build_features as wp22_build_features # noqa: E402222324def month_add(m: str, k: int) -> str:25 y, mm = int(m[:4]), int(m[5:7]); mm += k26 while mm > 12:27 mm -= 12; y += 128 while mm < 1:29 mm += 12; y -= 130 return f"{y:04d}-{mm:02d}"313233def months_between(a: str, b: str) -> int:34 return (int(b[:4]) - int(a[:4])) * 12 + int(b[5:7]) - int(a[5:7])353637def load_sales() -> pl.DataFrame:38 """WP22 analytic layer + WP22 feature engineering (imputed logs, age bins, grid cells)."""39 a = pl.read_parquet(C.WP22_ANALYTIC)40 f = wp22_build_features(a)41 keep = ["tx_id", "parcel_id", "date", "month", "year", "price", "log_price", "property_class", "subtype_c", "genre_c", "storeys_c", "units_c",42 "lfa", "lfa2", "lla", "age_f", "fa_missing", "la_missing", "age_missing", "fa_f", "la_f", "x_m", "y_m", "mun_code", "mrc_code", "region_code",43 "grid1k", "grid2k", "grid5k", "assessment_at_sale", "new_construction", "mamh_population"]44 return f.select([c for c in keep if c in f.columns]).sort("date")454647# ---- market features at an origin (real-time, from sales ≤ origin) ---------------------------------48def market_features_at(sales: pl.DataFrame, origin: str, key: str = "mun_code", prefix: str = "mun") -> pl.DataFrame:49 """Per key: log median price per m² level and momentum (3/6/12 m), volumes, dispersion, computed from sales ≤ origin.50 Price per m² is a crude real-time proxy of the local level; the quality-adjusted alternative (vintage index) is added in Task B."""51 s = sales.filter(pl.col("month") <= origin).with_columns((pl.col("log_price") - pl.col("lfa") * 0.5).alias("_q")) # crude quality strip52 out = None53 for lab, k in (("3", 3), ("6", 6), ("12", 12), ("24", 24)):54 w = s.filter(pl.col("month") > month_add(origin, -k))55 g = w.group_by(key).agg(pl.col("log_price").median().alias(f"{prefix}_medlp_{lab}"), pl.col("_q").median().alias(f"{prefix}_q_{lab}"),56 pl.len().alias(f"{prefix}_n_{lab}"), pl.col("log_price").std().alias(f"{prefix}_sd_{lab}"))57 out = g if out is None else out.join(g, on=key, how="full", coalesce=True)58 out = out.with_columns(59 (pl.col(f"{prefix}_q_3") - pl.col(f"{prefix}_q_12")).alias(f"{prefix}_mom_3_12"),60 (pl.col(f"{prefix}_q_6") - pl.col(f"{prefix}_q_24")).alias(f"{prefix}_mom_6_24"),61 (pl.col(f"{prefix}_n_3") / (pl.col(f"{prefix}_n_12") / 4)).alias(f"{prefix}_vol_ratio"),62 )63 return out646566def knn_features_at(sales: pl.DataFrame, origin: str, targets: pl.DataFrame, k: int = 10, window_months: int = 12) -> pl.DataFrame:67 """For each target sale: mean quality-stripped log price and distance of the k nearest sales of the same class dated in68 (origin-12 m, origin]. Targets may be dated after the origin; only past sales are used."""69 past = sales.filter((pl.col("month") <= origin) & (pl.col("month") > month_add(origin, -window_months)))70 rows = []71 for cls in C.PROPERTY_TYPES:72 p = past.filter(pl.col("property_class") == cls); t = targets.filter(pl.col("property_class") == cls)73 if p.height < k or t.height == 0:74 continue75 tree = cKDTree(p.select(["x_m", "y_m"]).to_numpy())76 q = (p["log_price"] - 0.5 * p["lfa"]).to_numpy()77 d, idx = tree.query(t.select(["x_m", "y_m"]).to_numpy(), k=k)78 rows.append(pl.DataFrame({"tx_id": t["tx_id"], "knn_q": q[idx].mean(axis=1), "knn_dist_km": d.mean(axis=1) / 1000, "knn_q_sd": q[idx].std(axis=1)}))79 schema = {"tx_id": pl.Utf8, "knn_q": pl.Float64, "knn_dist_km": pl.Float64, "knn_q_sd": pl.Float64}80 return pl.concat([r.cast(schema) for r in rows]) if rows else pl.DataFrame(schema=schema)818283def task_a_sets(sales: pl.DataFrame, origin: str, horizons=C.HORIZONS, train_months: int | None = None):84 """Returns (train, test) frames with features at the origin. train_months=None → expanding window."""85 tr = sales.filter(pl.col("month") <= origin)86 if train_months:87 tr = tr.filter(pl.col("month") > month_add(origin, -train_months))88 te = sales.filter((pl.col("month") > origin) & (pl.col("month") <= month_add(origin, max(horizons)))).with_columns(89 pl.col("month").map_elements(lambda m: months_between(origin, m), return_dtype=pl.Int32).alias("h"))90 mf = market_features_at(sales, origin, "mun_code", "mun"); rf = market_features_at(sales, origin, "region_code", "reg")91 # training rows get market features at *their own* month-end? No: to mimic forecasting, training features are computed92 # at each training sale's month (information ≤ that month), test features at the origin.93 def attach(df, at_origin: bool):94 if at_origin:95 d = df.join(mf, on="mun_code", how="left").join(rf, on="region_code", how="left")96 kn = knn_features_at(sales, origin, df)97 return d.join(kn, on="tx_id", how="left").with_columns(pl.lit(months_between("2021-01", origin)).alias("origin_t"))98 # training rows: each sale is assigned a horizon h ∈ HORIZONS (deterministic hash) and receives the market and99 # neighbour features of the origin h months before its sale month, so that the models learn horizon-dependence100 rng = np.random.default_rng(C.SEED); hs = np.array(horizons)101 df = df.with_columns(pl.Series("h", hs[rng.integers(0, len(hs), df.height)]).cast(pl.Int32))102 parts = []103 for (m, h), sub in df.group_by(["month", "h"]):104 mo = month_add(m, -int(h))105 if mo < "2021-01":106 continue107 d = sub.join(market_features_at(sales, mo, "mun_code", "mun"), on="mun_code", how="left").join(market_features_at(sales, mo, "region_code", "reg"), on="region_code", how="left")108 d = d.join(knn_features_at(sales, mo, sub), on="tx_id", how="left").with_columns(pl.lit(months_between("2021-01", mo)).alias("origin_t"))109 parts.append(d)110 return pl.concat(parts)111 return attach(tr, False), attach(te, True)112113114FEATURES_HEDONIC = ["lfa", "lfa2", "lla", "age_f", "fa_missing", "la_missing", "age_missing", "x_m", "y_m"]115FEATURES_CAT = ["property_class", "subtype_c", "genre_c", "storeys_c", "units_c", "region_code"]116FEATURES_MARKET = ["mun_medlp_3", "mun_q_3", "mun_q_12", "mun_n_3", "mun_n_12", "mun_sd_12", "mun_mom_3_12", "mun_mom_6_24", "mun_vol_ratio",117 "reg_q_3", "reg_q_12", "reg_mom_3_12", "reg_mom_6_24", "knn_q", "knn_dist_km", "knn_q_sd", "origin_t"]118119120def design(df: pl.DataFrame, cats: dict | None = None):121 """Numeric design matrix (one-hot categoricals with a fixed category list) + the target vector."""122 d = df.with_columns([pl.col(c).fill_null(0.0) for c in FEATURES_MARKET if c in df.columns])123 X = d.select([c for c in FEATURES_HEDONIC + FEATURES_MARKET if c in d.columns]).to_numpy().astype(np.float64)124 cols = [c for c in FEATURES_HEDONIC + FEATURES_MARKET if c in d.columns]125 cats = cats or {c: sorted(d[c].cast(pl.Utf8).fill_null("na").unique().to_list()) for c in FEATURES_CAT}126 oh = []127 for c, levels in cats.items():128 v = d[c].cast(pl.Utf8).fill_null("na").to_numpy()129 for lv in levels:130 oh.append((v == lv).astype(np.float64)); cols.append(f"{c}={lv}")131 X = np.column_stack([X] + oh) if oh else X132 X = np.nan_to_num(X)133 return X, d["log_price"].to_numpy(), cols, cats134