1"""Replication entry point for F6-conf-02. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3Reads only data/, writes out/results.json, out/figure1.png, out/table1.csv, out/exclusions.csv.4No randomness is used (random_seed = null)."""5import json, os6import numpy as np7import pandas as pd8from scipy import stats9import matplotlib10matplotlib.use("Agg")11import matplotlib.pyplot as plt1213DATA = "data/gc5y_daily.csv"14OUT = "out"15PRE_N = 12016CUT_DATE = "2020-03-04"17os.makedirs(OUT, exist_ok=True)1819raw = pd.read_csv(DATA)20n_raw = len(raw)21d = raw.copy()22d["date_dt"] = pd.to_datetime(d["date"], format="%Y-%m-%d")23issues = []24excl = []2526# ---- audit ----27n_na = int(d[["date", "yield_5y", "post_cut"]].isna().any(axis=1).sum())28n_dup = int(d["date"].duplicated().sum())29monotone = bool(d["date_dt"].is_monotonic_increasing)30n_treated = int((d["post_cut"] == 1).sum())31bad_flag = int((~d["post_cut"].isin([0, 1])).sum())32flag_mismatch = int(((d["date_dt"] >= CUT_DATE).astype(int) != d["post_cut"]).sum())33gaps = d["date_dt"].diff().dt.days34n_gaps_gt3 = int((gaps > 3).sum())3536if n_treated == 0:37 issues.append({38 "issue": f"Coverage: file spans {d['date'].iloc[0]} to {d['date'].iloc[-1]} (retrieval end_date=2019-12-31 per data/README.md); "39 f"it ends before the 2020-03-04 cut and contains 0 rows with post_cut = 1.",40 "handling": "Treated window is empty, so the post-cut mean, the pre-window (120 days before the first flagged day) and the effect "41 "are undefined. Declared premise_invalid / infeasible; inestimable keys set to null. No proxy (e.g. placebo date) substituted."})42issues.append({"issue": f"Missing values: {n_na} rows with missing fields; duplicate dates: {n_dup}; dates monotone increasing: {monotone}.",43 "handling": "No action needed; no rows removed."})44issues.append({"issue": f"post_cut values outside {{0,1}}: {bad_flag}; rows where post_cut disagrees with date >= {CUT_DATE}: {flag_mismatch}.",45 "handling": "post_cut used exactly as provided (not redefined)."})46issues.append({"issue": f"Calendar gaps: {n_gaps_gt3} consecutive-row gaps longer than 3 calendar days (max {int(gaps.max())} days), consistent with holidays/non-publication days.",47 "handling": "Per task rule, one-day changes are taken between consecutive rows regardless of calendar gap."})4849# ---- changes and baseline ----50d["chg_bp"] = 100.0 * d["yield_5y"].diff()51chg = d["chg_bp"].dropna()52excl.append(dict(step=1, rule="first row has no previous row: one-day change undefined", n_before=n_raw,53 n_after=len(chg), n_removed=n_raw - len(chg), rationale="change = 100*(y_t - y_{t-1}) needs a prior row"))54baseline_metric = float(chg.mean())5556# ---- effect (only if treated rows exist) ----57effect = lo = hi = model_metric = None58n_analysed = None59if n_treated > 0:60 first = int(np.flatnonzero(d["post_cut"].values == 1)[0])61 post = d.loc[d["post_cut"] == 1, "chg_bp"].dropna()62 pre = d["chg_bp"].iloc[max(first - PRE_N, 0):first].dropna()63 if len(pre) >= 2 and len(post) >= 2:64 effect = float(post.mean() - pre.mean())65 v1, v0 = post.var(ddof=1) / len(post), pre.var(ddof=1) / len(pre)66 se = np.sqrt(v1 + v0)67 df = (v1 + v0) ** 2 / (v1 ** 2 / (len(post) - 1) + v0 ** 2 / (len(pre) - 1))68 t = stats.t.ppf(0.975, df)69 lo, hi = float(effect - t * se), float(effect + t * se)70 resid = np.concatenate([post - post.mean(), pre - pre.mean()])71 model_metric = float(np.sqrt(np.mean(resid ** 2)))72 n_analysed = int(len(post) + len(pre))73 excl.append(dict(step=2, rule="keep 120 pre rows before first flagged day + all post_cut=1 rows", n_before=len(chg),74 n_after=n_analysed, n_removed=len(chg) - n_analysed, rationale="estimand windows"))75else:76 excl.append(dict(step=2, rule="keep 120 pre rows before first flagged day + all post_cut=1 rows", n_before=len(chg),77 n_after=0, n_removed=len(chg), rationale="no post_cut=1 rows: windows undefined, nothing analysed"))7879if effect is None:80 conclusion, feas = "premise_invalid", "infeasible"81 reason = ("data/gc5y_daily.csv ends on " + d["date"].iloc[-1] + " and has no rows with post_cut = 1; the post-cut window "82 "and the 120-day pre-window defined relative to the first flagged day do not exist, so the effect cannot be computed.")83else:84 conclusion = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")85 feas, reason = "complete", None8687res = {88 "task_id": "F6-conf-02",89 "sample_start": d["date"].iloc[0], "sample_end": d["date"].iloc[-1],90 "n_treated": n_treated,91 "baseline_metric": baseline_metric,92 "model_metric": model_metric,93 "effect_estimate": effect, "effect_ci_low": lo, "effect_ci_high": hi,94 "conclusion": conclusion, "feasibility": feas, "feasibility_reason": reason,95 "data_issues": issues,96 "n_input_rows": {"gc5y_daily.csv": n_raw},97 "n_analysed": n_analysed,98 "random_seed": None,99 "baseline_n_changes": int(len(chg)),100 "baseline_sd_change_bp": float(chg.std(ddof=1)),101}102with open(f"{OUT}/results.json", "w") as f:103 json.dump(res, f, indent=2)104pd.DataFrame(excl, columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv(f"{OUT}/exclusions.csv", index=False)105106# ---- table ----107d["year"] = d["date_dt"].dt.year108tab = d.groupby("year").agg(n_days=("yield_5y", "size"), mean_yield=("yield_5y", "mean"),109 mean_change_bp=("chg_bp", "mean")).reset_index()110tab.round(4).to_csv(f"{OUT}/table1.csv", index=False)111112# ---- figure ----113fig, ax = plt.subplots(figsize=(9, 4))114ax.plot(d["date_dt"], d["yield_5y"], lw=0.8)115ax.set_title("Government of Canada 5-year benchmark yield (BD.CDN.5YR.DQ.YLD)\nFile ends "116 + d["date"].iloc[-1] + "; no post_cut = 1 rows (cut date 2020-03-04 not covered)", fontsize=9)117ax.set_ylabel("Yield (%)"); ax.set_xlabel("Date"); ax.grid(alpha=0.3)118fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120)119120print(f"sample {res['sample_start']}..{res['sample_end']}, n_treated={n_treated}, baseline={baseline_metric:.4f} bp, "121 f"conclusion={conclusion}, feasibility={feas}")122