1"""Replication entry point for F6-conf-02. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3Reads only data/gc5y_daily.csv; writes out/results.json, out/table1.csv,4out/figure1.png, out/exclusions.csv. No randomness is used."""5import json, os6import numpy as np7import pandas as pd8from scipy import stats9import matplotlib10matplotlib.use("Agg")11import matplotlib.pyplot as plt1213TASK_ID = "F6-conf-02"14FN = "gc5y_daily.csv"15PRE_N = 12016os.makedirs("out", exist_ok=True)1718raw = pd.read_csv(os.path.join("data", FN))19n_raw = len(raw)20df = raw.copy()21df["date"] = pd.to_datetime(df["date"], format="%Y-%m-%d")22df = df.sort_values("date").reset_index(drop=True)2324# ---- audit ----25issues = []26n_missing = int(df.isna().sum().sum())27n_dup = int(df["date"].duplicated().sum())28if n_missing:29 issues.append({"issue": f"{n_missing} missing values", "handling": "rows with missing yield dropped"})30if n_dup:31 issues.append({"issue": f"{n_dup} duplicate dates", "handling": "reported; not removed"})32n_treated = int((df["post_cut"] == 1).sum())33sample_start = df["date"].min().strftime("%Y-%m-%d")34sample_end = df["date"].max().strftime("%Y-%m-%d")35if n_treated == 0:36 issues.append({37 "issue": f"No rows with post_cut = 1: the file covers {sample_start} to {sample_end} "38 "(README retrieval used end_date=2019-12-31), entirely before the 2020-03-04 cut; "39 "post_cut is 0 for all rows.",40 "handling": "Estimand not computable (no treated window, and the 120-day pre-window is defined "41 "relative to the first flagged day, which does not exist). post_cut not redefined; "42 "no proxy estimate substituted; effect keys set to null, conclusion premise_invalid."})43mismatch = int(((df["date"] >= "2020-03-04") != (df["post_cut"] == 1)).sum())44if mismatch:45 issues.append({"issue": f"{mismatch} rows where post_cut disagrees with its documented definition",46 "handling": "post_cut used exactly as provided, per task rules"})4748# ---- changes and baseline ----49df["chg_bp"] = 100.0 * df["yield_5y"].diff()50chg = df["chg_bp"].dropna()51baseline_metric = float(chg.mean())52excl = [{"step": 1, "rule": "first row has no previous row: one-day change undefined",53 "n_before": n_raw, "n_after": len(chg), "n_removed": n_raw - len(chg),54 "rationale": "changes are between consecutive rows (TASK.md)"}]5556# ---- main estimate ----57effect = lo = hi = rmse = None58n_analysed = None59feas, feas_reason, conclusion = "complete", None, None60if n_treated > 0:61 first = df.index[df["post_cut"] == 1][0]62 post = df.loc[df["post_cut"] == 1, "chg_bp"].dropna()63 pre = df.loc[max(first - PRE_N, 0):first - 1, "chg_bp"].dropna()64 if len(pre) >= 2 and len(post) >= 2:65 effect = float(post.mean() - pre.mean())66 v1, v0 = post.var(ddof=1) / len(post), pre.var(ddof=1) / len(pre)67 se = np.sqrt(v1 + v0)68 dfw = (v1 + v0) ** 2 / (v1 ** 2 / (len(post) - 1) + v0 ** 2 / (len(pre) - 1))69 t = stats.t.ppf(0.975, dfw)70 lo, hi = float(effect - t * se), float(effect + t * se)71 resid = np.concatenate([post - post.mean(), pre - pre.mean()])72 rmse = float(np.sqrt(np.mean(resid ** 2)))73 n_analysed = int(len(post) + len(pre))74 conclusion = "positive_effect" if lo > 0 else "negative_effect" if hi < 0 else "no_detectable_effect"75 excl.append({"step": 2, "rule": "keep 120 pre rows + post_cut=1 rows",76 "n_before": len(chg), "n_after": n_analysed, "n_removed": len(chg) - n_analysed,77 "rationale": "estimand windows"})78if conclusion is None:79 feas, conclusion = "infeasible", "premise_invalid"80 feas_reason = (f"The file contains no rows with post_cut = 1 (sample {sample_start} to {sample_end}, "81 "ending before the 2020-03-04 cut), so the post-cut mean change and the 120-day "82 "pre-window anchored at the first flagged day cannot be computed.")83 excl.append({"step": 2, "rule": "keep 120 rows before first post_cut=1 row + post_cut=1 rows",84 "n_before": len(chg), "n_after": 0, "n_removed": len(chg),85 "rationale": "no post_cut=1 rows exist, so both estimand windows are empty"})86pd.DataFrame(excl).to_csv("out/exclusions.csv", index=False)8788# ---- table and figure ----89df["year"] = df["date"].dt.year90tab = df.groupby("year").agg(n_days=("yield_5y", "size"), mean_yield=("yield_5y", "mean"),91 mean_change_bp=("chg_bp", "mean")).reset_index()92tab.round(6).to_csv("out/table1.csv", index=False)9394fig, ax = plt.subplots(figsize=(9, 4))95ax.plot(df["date"], df["yield_5y"], lw=0.8, color="tab:blue")96ax.axvline(pd.Timestamp("2020-03-04"), color="red", ls="--", lw=1, label="2020-03-04 cut (outside sample)")97ax.set_xlim(df["date"].min(), pd.Timestamp("2020-06-30"))98ax.set_ylabel("5-year GoC yield (%)")99ax.set_title(f"GoC 5-year benchmark yield, {sample_start} to {sample_end} (no post_cut=1 rows)"100 if n_treated == 0 else "GoC 5-year benchmark yield")101ax.legend(loc="upper right")102fig.tight_layout()103fig.savefig("out/figure1.png", dpi=120)104105res = {106 "task_id": TASK_ID, "sample_start": sample_start, "sample_end": sample_end,107 "n_treated": n_treated, "baseline_metric": baseline_metric, "model_metric": rmse,108 "effect_estimate": effect, "effect_ci_low": lo, "effect_ci_high": hi,109 "conclusion": conclusion, "feasibility": feas, "feasibility_reason": feas_reason,110 "data_issues": issues, "n_input_rows": {FN: n_raw}, "n_analysed": n_analysed,111 "random_seed": None,112 "baseline_n_changes": int(len(chg)),113}114with open("out/results.json", "w") as f:115 json.dump(res, f, indent=2)116print(f"rows={n_raw} treated={n_treated} baseline_mean_change_bp={baseline_metric:.6f} conclusion={conclusion}")117