1"""F6-conf-02 replication script. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3Regenerates out/results.json, out/figure1.png, out/table1.csv from data/ only."""4import json, os5import numpy as np, pandas as pd6import matplotlib7matplotlib.use("Agg")8import matplotlib.pyplot as plt910SEED = None # no randomness used11os.makedirs("out", exist_ok=True)12d = pd.read_csv("data/gc5y_daily.csv")13n_rows = len(d)14d["date_dt"] = pd.to_datetime(d["date"], format="%Y-%m-%d")1516# ---- audit ----17issues = []18n_missing = int(d[["date", "yield_5y", "post_cut"]].isna().sum().sum())19n_dup = int(d["date"].duplicated().sum())20mono = bool(d["date_dt"].is_monotonic_increasing)21gaps = d["date_dt"].diff().dt.days22n_treated = int((d["post_cut"] == 1).sum())23expected = (d["date_dt"] >= "2020-03-04").astype(int)24n_flag_mismatch = int((expected != d["post_cut"]).sum())25max_date = d["date"].max()2627if n_treated == 0:28 issues.append({29 "issue": f"Sample ends {max_date} (README retrieval end_date=2019-12-31), before the 2020-03-04 cut; "30 "post_cut = 0 for all rows, so there are zero treated days and no post-cut window.",31 "handling": "Estimand not computable; effect, CI and model RMSE set to null; conclusion premise_invalid. "32 "No proxy (e.g. other dates or synthetic data) substituted."})33issues.append({34 "issue": "README describes 2020-03-04 both as the announcement date and as 'the first business day after' the announcement (internal inconsistency in flag description).",35 "handling": "Moot because no 2020 rows exist; post_cut used exactly as provided, not redefined."})36issues.append({37 "issue": f"Missing values: {n_missing}; duplicate dates: {n_dup}; dates monotonic increasing: {mono}; "38 f"post_cut disagreements with 'date >= 2020-03-04' rule: {n_flag_mismatch}.",39 "handling": "No cleaning required; recorded for transparency."})40issues.append({41 "issue": f"Calendar gaps between consecutive rows: max {int(gaps.max())} days; "42 f"{int((gaps > 3).sum())} gaps longer than a weekend (Canadian holidays / non-publication days).",43 "handling": "Per task rules, one-day changes are taken between consecutive rows regardless of calendar gap."})4445# ---- descriptive baseline ----46d["chg_bp"] = 100 * d["yield_5y"].diff()47baseline = float(d["chg_bp"].mean())4849# ---- main estimate (only if treated rows exist) ----50effect = lo = hi = rmse = None51n_analysed = None52if n_treated > 0:53 from scipy import stats54 first = d.index[d["post_cut"] == 1][0]55 post = d.loc[d["post_cut"] == 1, "chg_bp"].dropna()56 pre = d.loc[max(first - 120, 1):first - 1, "chg_bp"].dropna()57 if len(pre) == 120 and len(post) >= 2:58 effect = float(post.mean() - pre.mean())59 v1, v0 = post.var(ddof=1) / len(post), pre.var(ddof=1) / len(pre)60 se = np.sqrt(v1 + v0)61 df = (v1 + v0) ** 2 / (v1**2 / (len(post) - 1) + v0**2 / (len(pre) - 1))62 t = stats.t.ppf(0.975, df)63 lo, hi = effect - t * se, effect + t * se64 resid = np.concatenate([post - post.mean(), pre - pre.mean()])65 rmse = float(np.sqrt(np.mean(resid**2)))66 n_analysed = int(len(post) + len(pre))6768if effect is None:69 conclusion, feas = "premise_invalid", "infeasible"70 reason = (f"The file covers {d['date'].min()} to {max_date} only; n_treated = 0, so there is no post-cut "71 "window and the post-minus-pre difference in mean one-day changes cannot be computed.")72else:73 conclusion = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")74 feas, reason = "complete", None7576res = {77 "task_id": "F6-conf-02",78 "sample_start": d["date"].min(), "sample_end": max_date,79 "n_treated": n_treated,80 "baseline_metric": round(baseline, 6),81 "model_metric": rmse, "effect_estimate": effect,82 "effect_ci_low": lo, "effect_ci_high": hi,83 "conclusion": conclusion, "feasibility": feas, "feasibility_reason": reason,84 "data_issues": issues,85 "n_input_rows": {"gc5y_daily.csv": n_rows},86 "n_analysed": n_analysed, "random_seed": SEED,87}88with open("out/results.json", "w") as f:89 json.dump(res, f, indent=2)9091# ---- table ----92d["year"] = d["date_dt"].dt.year93tab = d.groupby("year").agg(n_days=("yield_5y", "size"), mean_yield=("yield_5y", "mean"),94 mean_change_bp=("chg_bp", "mean")).reset_index()95tab.round(4).to_csv("out/table1.csv", index=False)9697# ---- figure ----98fig, ax = plt.subplots(figsize=(9, 4))99ax.plot(d["date_dt"], d["yield_5y"], lw=0.8)100ax.set_title(f"GoC 5-year benchmark yield, {d['date'].min()} to {max_date}\n"101 "(no observations on/after 2020-03-04 in the file)")102ax.set_ylabel("Yield (%)"); ax.set_xlabel("Date")103fig.tight_layout(); fig.savefig("out/figure1.png", dpi=120)104print(json.dumps(res, indent=2))105