1"""Replication entry point for F6-conf-02. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3"""4import json, os5import numpy as np, pandas as pd6import matplotlib7matplotlib.use("Agg")8import matplotlib.pyplot as plt9from scipy import stats1011SEED = None # no randomness used12os.makedirs("out", exist_ok=True)13raw = pd.read_csv("data/gc5y_daily.csv")14n_rows = len(raw)15d = raw.copy()16d["date"] = pd.to_datetime(d["date"], errors="coerce")17issues = []1819# ---- audit ----20if d["date"].isna().any():21 issues.append({"issue": f"{int(d['date'].isna().sum())} unparseable dates", "handling": "dropped"})22if d["yield_5y"].isna().any():23 issues.append({"issue": f"{int(d['yield_5y'].isna().sum())} missing yields", "handling": "dropped before differencing"})24if d["date"].duplicated().any():25 issues.append({"issue": "duplicated dates", "handling": "kept first"})26d = d.dropna(subset=["date", "yield_5y"]).drop_duplicates("date").sort_values("date").reset_index(drop=True)27d["chg_bp"] = 100 * d["yield_5y"].diff()2829n_treated = int((d["post_cut"] == 1).sum())30last = d["date"].max()31if n_treated == 0:32 issues.append({33 "issue": (f"File covers {d['date'].min().date()} to {last.date()} only (README retrieval end_date=2019-12-31); "34 "zero rows have post_cut=1 and no observation exists on/after the 2020-03-04 cut, "35 "so the post-cut window is empty and the 120-day pre-window (defined relative to the first flagged day) is undefined."),36 "handling": "Estimand declared not estimable (premise_invalid); effect, interval, model RMSE and n_analysed set to null. No proxy window substituted."})37mism = int(((d["date"] >= "2020-03-04").astype(int) != d["post_cut"]).sum())38if mism:39 issues.append({"issue": f"{mism} rows where post_cut disagrees with date>=2020-03-04",40 "handling": "post_cut used exactly as provided per task rules"})41gaps = d["date"].diff().dt.days42ng = int((gaps > 4).sum())43if ng:44 issues.append({"issue": f"{ng} gaps of >4 calendar days between consecutive rows (holidays/non-publication days)",45 "handling": "changes computed between consecutive rows as specified; multi-day gaps retained"})4647baseline = float(d["chg_bp"].mean())4849effect = lo = hi = rmse = n_an = None50if n_treated > 0:51 first = d.index[d["post_cut"] == 1][0]52 post = d.loc[d["post_cut"] == 1, "chg_bp"].dropna()53 pre = d.loc[max(first - 120, 1):first - 1, "chg_bp"].dropna()54 if len(pre) >= 2 and len(post) >= 2:55 effect = float(post.mean() - pre.mean())56 v1, v0 = post.var(ddof=1) / len(post), pre.var(ddof=1) / len(pre)57 se = np.sqrt(v1 + v0)58 df = (v1 + v0) ** 2 / (v1**2 / (len(post) - 1) + v0**2 / (len(pre) - 1))59 t = stats.t.ppf(0.975, df)60 lo, hi = effect - t * se, effect + t * se61 resid = np.concatenate([post - post.mean(), pre - pre.mean()])62 rmse = float(np.sqrt(np.mean(resid**2)))63 n_an = int(len(post) + len(pre))6465if effect is None:66 concl, feas = "premise_invalid", "infeasible"67 reason = ("No rows with post_cut=1: the file ends on 2019-12-31, before the 2020-03-04 cut, "68 "so neither the post-cut mean nor the 120-day pre-window can be formed.")69else:70 concl = "positive_effect" if lo > 0 else "negative_effect" if hi < 0 else "no_detectable_effect"71 feas, reason = "complete", None7273res = {"task_id": "F6-conf-02",74 "sample_start": str(d["date"].min().date()), "sample_end": str(last.date()),75 "n_treated": n_treated, "baseline_metric": baseline, "model_metric": rmse,76 "effect_estimate": effect, "effect_ci_low": None if lo is None else float(lo),77 "effect_ci_high": None if hi is None else float(hi),78 "conclusion": concl, "feasibility": feas, "feasibility_reason": reason,79 "data_issues": issues, "n_input_rows": {"gc5y_daily.csv": n_rows},80 "n_analysed": n_an, "random_seed": SEED}81json.dump(res, open("out/results.json", "w"), indent=2)8283# table84d["year"] = d["date"].dt.year85tab = d.groupby("year").agg(n_days=("yield_5y", "size"), mean_yield=("yield_5y", "mean"),86 mean_change_bp=("chg_bp", "mean")).reset_index()87tab.round(4).to_csv("out/table1.csv", index=False)8889# figure90fig, ax = plt.subplots(figsize=(9, 4))91ax.plot(d["date"], d["yield_5y"], lw=0.8)92ax.axvline(pd.Timestamp("2020-03-04"), color="r", ls="--", lw=1, label="2020-03-04 cut (outside sample)")93ax.set_xlim(d["date"].min(), pd.Timestamp("2020-06-30"))94ax.set_ylabel("GoC 5-year yield (%)"); ax.set_title("GoC 5-year benchmark yield, file sample (ends 2019-12-31)")95ax.legend(); fig.tight_layout(); fig.savefig("out/figure1.png", dpi=120)96print(json.dumps(res, indent=2))97