1"""H3 — Concentration and the long tail of merchant catalogue sizes (S_A, all stores with >= 1 listing).2Gini, top-1 % / top-10 % shares (store bootstrap CIs), rank-size slope, discrete power-law vs log-normal (CSN MLE + Vuong LR).3Within-region Gini (S_R). Preregistered decision: Gini >= 0.70 and top-10 % share >= 0.50.4"""5from __future__ import annotations6import json7import numpy as np, pandas as pd8from scipy import stats, optimize, special9from src.common import ROOT, SEED, get_logger, ledger, write_json, PROC_DIR, EXP_DIR, P1011log = get_logger("models.h3")12PROC = PROC_DIR131415def gini(x):16 x = np.sort(np.asarray(x, float)); n = len(x)17 return float((2 * np.arange(1, n + 1) - n - 1).dot(x) / (n * x.sum()))181920def top_share(x, q):21 x = np.sort(np.asarray(x, float))[::-1]; k = max(1, int(np.ceil(q * len(x)))); return float(x[:k].sum() / x.sum())222324def powerlaw_fit(x):25 """Discrete power law with x_min chosen by KS minimisation (Clauset, Shalizi & Newman 2009)."""26 x = np.asarray(x, int); best = None27 for xmin in np.unique(x)[:-10]:28 t = x[x >= xmin]; n = len(t)29 if n < 50: break30 nll = lambda a: -(n * -np.log(special.zeta(a, xmin)) - a * np.log(t).sum()) if a > 1.0001 else 1e1831 r = optimize.minimize_scalar(nll, bounds=(1.01, 6), method="bounded"); a = r.x32 # KS distance33 xs = np.arange(xmin, t.max() + 1); cdf_m = np.cumsum(xs ** -a) / special.zeta(a, xmin)34 emp = np.searchsorted(np.sort(t), xs, side="right") / n35 ks = np.max(np.abs(emp - cdf_m))36 if best is None or ks < best["ks"]:37 best = {"xmin": int(xmin), "alpha": float(a), "ks": float(ks), "n_tail": int(n), "ll_pl": float(-nll(a))}38 return best394041def lognormal_vs_powerlaw(x, fit):42 t = np.asarray(x, float); t = t[t >= fit["xmin"]]43 # discrete log-normal truncated at xmin: approximate by continuous log-normal density ratio on the same tail44 mu, sig = np.log(t).mean(), np.log(t).std(ddof=0)45 ll_ln = stats.lognorm.logpdf(t, s=sig, scale=np.exp(mu)) - np.log(1 - stats.lognorm.cdf(fit["xmin"] - 0.5, s=sig, scale=np.exp(mu)))46 a, xmin = fit["alpha"], fit["xmin"]47 ll_pl = -a * np.log(t) - np.log(special.zeta(a, xmin))48 d = ll_pl - ll_ln; R = d.sum(); se = d.std(ddof=1) * np.sqrt(len(d))49 return {"vuong_R": float(R), "vuong_z": float(R / se) if se > 0 else None, "vuong_p": float(2 * (1 - stats.norm.cdf(abs(R / se)))) if se > 0 else None,50 "favoured": "power-law" if R > 0 else "log-normal", "lognormal_mu": float(mu), "lognormal_sigma": float(sig)}515253def main():54 sp = pd.read_parquet(PROC / "store_panel.parquet"); rng = np.random.default_rng(SEED)55 x = sp.n_listings.values56 B = 2000; g_b = np.empty(B); t1 = np.empty(B); t10 = np.empty(B)57 for b in range(B):58 d = rng.choice(x, len(x), replace=True); g_b[b] = gini(d); t1[b] = top_share(d, 0.01); t10[b] = top_share(d, 0.10)59 ranks = np.arange(1, len(x) + 1); xs = np.sort(x)[::-1]60 slope = float(np.polyfit(np.log(ranks), np.log(xs), 1)[0])61 fit = powerlaw_fit(x); lv = lognormal_vs_powerlaw(x, fit)62 res = {"stores": int(len(x)), "listings": int(x.sum()), "gini": round(gini(x), 4), "gini_ci": [round(float(np.percentile(g_b, 2.5)), 4), round(float(np.percentile(g_b, 97.5)), 4)],63 "top1pct_share": round(top_share(x, 0.01), 4), "top1_ci": [round(float(np.percentile(t1, 2.5)), 4), round(float(np.percentile(t1, 97.5)), 4)],64 "top10pct_share": round(top_share(x, 0.10), 4), "top10_ci": [round(float(np.percentile(t10, 2.5)), 4), round(float(np.percentile(t10, 97.5)), 4)],65 "median_catalogue": float(np.median(x)), "p90": float(np.percentile(x, 90)), "p99": float(np.percentile(x, 99)), "max": int(x.max()),66 "rank_size_slope_loglog": round(slope, 4), "powerlaw": fit, "lognormal_vs_powerlaw": lv,67 "decision_prereg": bool(gini(x) >= 0.70 and top_share(x, 0.10) >= 0.50)}68 # within-region69 wr = {}70 for r, g in sp[sp.region_final != "Inconnue"].groupby("region_final"):71 v = g.n_listings.values72 if len(v) < 5: wr[r] = {"stores": int(len(v)), "gini": None}; continue73 gb = [gini(rng.choice(v, len(v), replace=True)) for _ in range(500)]74 wr[r] = {"stores": int(len(v)), "listings": int(v.sum()), "gini": round(gini(v), 3), "gini_ci": [round(float(np.percentile(gb, 2.5)), 3), round(float(np.percentile(gb, 97.5)), 3)],75 "top10pct_share": round(top_share(v, 0.10), 3), "median": float(np.median(v))}76 res["within_region"] = wr77 pd.DataFrame({"rank": ranks, "n_listings": xs}).to_csv(PROC / "rank_size.csv", index=False)78 write_json(res, EXP_DIR / "h3_longtail.json", "h3"); ledger("models.h3_longtail", **{k: v for k, v in res.items() if k != "within_region"})79 log.info(json.dumps({k: v for k, v in res.items() if k != "within_region"}))808182if __name__ == "__main__":83 main()84