1"""Listing-day risk set with time-varying local competition (Phase 3, analysis_plan §2).23Competitor set of listing i on day t: listings of the same make & model, model year4within ±1, at a *different* dealer, within radius R km, active on day t, drawn from5the reliable-source universe (S2L: direct dealers, full portals, banners; not6oscillating). Radii 10 / 25 / 50 km. Day index: UTC calendar days since the first7crawl day (2026-08-12 = day 0).8"""9from __future__ import annotations1011import json12import math13import time14from pathlib import Path1516import numpy as np17import polars as pl1819from . import config as C2021RADII = (10.0, 25.0, 50.0)22DAY0 = math.floor(C.WINDOW_START_TS / 86400)232425def day_of(ts: np.ndarray) -> np.ndarray:26 return (np.floor(ts / 86400) - DAY0).astype(np.int32)272829def haversine_km(lat1, lng1, lat2, lng2):30 R = 6371.031 p1, p2 = np.radians(lat1), np.radians(lat2)32 dphi = p2[None, :] - p1[:, None]33 dl = np.radians(lng2)[None, :] - np.radians(lng1)[:, None]34 a = np.sin(dphi / 2) ** 2 + np.cos(p1)[:, None] * np.cos(p2)[None, :] * np.sin(dl / 2) ** 235 return 2 * R * np.arcsin(np.sqrt(np.clip(a, 0, 1)))363738def daily_price_matrix(ev: pl.DataFrame, ids: list[str], e_day: np.ndarray, T: int) -> tuple[np.ndarray, np.ndarray]:39 """Step-function log price per listing-day (nan when no price yet) and daily cut indicator."""40 n = len(ids)41 P = np.full((n, T), np.nan, dtype=np.float32)42 K = np.zeros((n, T), dtype=np.float32)43 idx = {l: i for i, l in enumerate(ids)}44 g = ev.filter(pl.col("listing_id").is_in(ids)).sort(["listing_id", "ts"]).group_by("listing_id", maintain_order=True).agg([pl.col("ts"), pl.col("price"), pl.col("kind")])45 days = np.arange(T)46 for lid, ts, pr, kd in g.iter_rows():47 i = idx[lid]48 d = day_of(np.asarray(ts, dtype=np.float64))49 d = np.minimum(d, T - 1)50 # price at end of each day = last event with day <= t51 pos = np.searchsorted(d, days, side="right") - 152 valid = pos >= 053 P[i, valid] = np.log(np.asarray(pr, dtype=np.float64))[pos[valid]]54 for dd, k in zip(d, kd):55 if k == "cut" and 0 <= dd < T:56 K[i, dd] = 1.057 return P, K585960def build(listings_full: Path = C.INTERIM / "listings_full.parquet", events: Path = C.PROCESSED / "price_events.parquet",61 out: Path = C.PROCESSED / "risk_set.parquet") -> dict:62 t0 = time.time()63 L = pl.read_parquet(listings_full)64 E = pl.read_parquet(events)65 T = int(day_of(np.array([C.SNAPSHOT_TS]))[0]) + 166 U = L.filter(pl.col("s2l_dynamic_with_stock") & pl.col("lat").is_not_null()).select(67 ["listing_id", "dealer_id", "make", "model", "year", "lat", "lng", "entry_ts", "end_ts", "exit_confirmed", "s2_dynamic", "left_censored", "first_cut_ts"])68 U = U.with_columns([69 pl.Series("e_day", day_of(U["entry_ts"].to_numpy())),70 pl.Series("x_day", np.minimum(day_of(U["end_ts"].to_numpy()), T - 1)),71 ])72 parts = []73 n_groups = 074 for (mk, md), G in U.group_by(["make", "model"]):75 n_groups += 176 ids = G["listing_id"].to_list()77 n = len(ids)78 e = G["e_day"].to_numpy(); x = G["x_day"].to_numpy(); y = G["year"].to_numpy()79 lat = G["lat"].to_numpy(); lng = G["lng"].to_numpy()80 dealer = pl.Series(G["dealer_id"]).cast(pl.Categorical).to_physical().to_numpy()81 A = ((np.arange(T)[None, :] >= e[:, None]) & (np.arange(T)[None, :] <= x[:, None])).astype(np.float32)82 P, K = daily_price_matrix(E, ids, e, T)83 Pf = np.where(np.isnan(P), 0.0, P).astype(np.float32)84 has_p = (~np.isnan(P)).astype(np.float32) * A85 D = haversine_km(lat, lng, lat, lng)86 same_year = np.abs(y[:, None] - y[None, :]) <= 187 other_dealer = dealer[:, None] != dealer[None, :]88 base = same_year & other_dealer89 cols = {}90 for R in RADII:91 M = (base & (D <= R)).astype(np.float32)92 dens = M @ A # n x T competitor counts93 ccut = M @ K # competitor cuts on day t94 sp = M @ (Pf * has_p); cnt = M @ has_p95 mean_lp = np.where(cnt > 0, sp / np.maximum(cnt, 1), np.nan)96 tag = f"{int(R)}"97 cols[f"dens_{tag}"] = dens; cols[f"ccut_{tag}"] = ccut; cols[f"cmlp_{tag}"] = mean_lp98 # within-dealer same-model stock (own-dealer substitutes)99 Mo = (same_year & ~other_dealer).astype(np.float32) - np.eye(n, dtype=np.float32)100 own = Mo @ A101 # explode to listing-days102 rows_i, rows_t = [], []103 for i in range(n):104 tt = np.arange(e[i], x[i] + 1)105 rows_i.append(np.full(len(tt), i)); rows_t.append(tt)106 ii = np.concatenate(rows_i); tt = np.concatenate(rows_t)107 part = {"listing_id": np.array(ids, dtype=object)[ii], "t": tt.astype(np.int16), "listing_age": (tt - e[ii]).astype(np.int16),108 "logp": P[ii, tt], "cut": K[ii, tt].astype(np.int8), "own_stock": own[ii, tt].astype(np.int16)}109 for k, v in cols.items():110 part[k] = v[ii, tt]111 parts.append(pl.DataFrame(part))112 RS = pl.concat(parts)113 # exit indicator: last day of the spell for confirmed exits114 last = U.select(["listing_id", "x_day", "exit_confirmed", "e_day", "s2_dynamic", "left_censored", "first_cut_ts"])115 RS = RS.join(last, on="listing_id", how="left").with_columns([116 ((pl.col("t") == pl.col("x_day")) & pl.col("exit_confirmed")).cast(pl.Int8).alias("exit"),117 (pl.col("t") == pl.col("x_day")).cast(pl.Int8).alias("last_day"),118 ])119 # risk set for the first markdown: days up to and including the first cut day120 fc_day = pl.when(pl.col("first_cut_ts").is_not_null()).then(((pl.col("first_cut_ts") / 86400).floor() - DAY0).cast(pl.Int32)).otherwise(None)121 RS = RS.with_columns(fc_day.alias("first_cut_day"))122 RS = RS.with_columns([123 (pl.col("first_cut_day").is_null() | (pl.col("t") <= pl.col("first_cut_day"))).alias("in_first_cut_risk"),124 ((pl.col("t") == pl.col("first_cut_day"))).fill_null(False).cast(pl.Int8).alias("first_cut"), # NULL first_cut_day -> 0 (never cut), see E-016125 ])126 for tag in ("10", "25", "50"):127 RS = RS.with_columns((pl.col("logp") - pl.col(f"cmlp_{tag}")).alias(f"pricepos_{tag}"))128 RS = RS.drop(["x_day", "first_cut_ts"]).sort(["listing_id", "t"])129 RS.write_parquet(out, compression="zstd")130 stats = {"listing_days": RS.height, "listings": RS["listing_id"].n_unique(), "groups": n_groups, "T": T,131 "exits": int(RS["exit"].sum()), "first_cuts": int(RS["first_cut"].sum()), "seconds": round(time.time() - t0, 1)}132 (C.PROCESSED / "risk_set_stats.json").write_text(json.dumps(stats, indent=2))133 return stats134135136if __name__ == "__main__":137 print(json.dumps(build(), indent=2))138