1"""F6-conf-02 replication entry point. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3Reads only data/gc5y_daily.csv; 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 np, pandas as pd7from scipy import stats8import matplotlib9matplotlib.use("Agg")10import matplotlib.pyplot as plt1112DATA = "data/gc5y_daily.csv"13OUT = "out"14PRE_N = 12015os.makedirs(OUT, exist_ok=True)1617df = pd.read_csv(DATA, dtype={"date": str})18n_raw = len(df)19df["dt"] = pd.to_datetime(df["date"], format="%Y-%m-%d")20issues = []2122# ---- audit ----23n_missing = int(df[["date", "yield_5y", "post_cut"]].isna().sum().sum())24n_dup = int(df["date"].duplicated().sum())25if n_missing:26 issues.append({"issue": f"{n_missing} missing values", "handling": "reported; rows kept"})27if n_dup:28 issues.append({"issue": f"{n_dup} duplicate dates", "handling": "reported"})29if not df["dt"].is_monotonic_increasing:30 issues.append({"issue": "dates not sorted", "handling": "file order used as required by the change definition"})31bad_flag = int((~df["post_cut"].isin([0, 1])).sum())32if bad_flag:33 issues.append({"issue": f"{bad_flag} post_cut values not in {{0,1}}", "handling": "reported"})34n_treated = int((df["post_cut"] == 1).sum())35start, end = df["date"].iloc[0], df["date"].iloc[-1]36gaps = df["dt"].diff().dt.days37if gaps.max() > 5:38 issues.append({"issue": f"max calendar gap {int(gaps.max())} days", "handling": "reported"})3940# ---- one-day changes (bp), between consecutive rows, dated t ----41df["chg_bp"] = 100.0 * df["yield_5y"].diff()42chg = df.dropna(subset=["chg_bp"])43excl = [dict(step=1, rule="load file", n_before=n_raw, n_after=n_raw, n_removed=0,44 rationale="all rows read; no missing values or duplicates"),45 dict(step=2, rule="first difference (first row has no t-1)", n_before=n_raw,46 n_after=len(chg), n_removed=n_raw - len(chg),47 rationale="one-day change undefined for first row; remaining changes form baseline sample")]4849# ---- baseline: mean one-day change over whole file ----50baseline = float(chg["chg_bp"].mean())5152# ---- effect: post (post_cut=1) vs 120 trading days before first flagged day ----53effect = lo = hi = rmse = None54n_analysed = None55feas, reason, conclusion = "complete", None, None56if n_treated == 0:57 issues.append({"issue": f"no rows with post_cut = 1: the file ends {end}, before the 2020-03-04 cut "58 "(README: retrieved with end_date=2019-12-31)",59 "handling": "treated window empty and pre-window undefined; estimand declared not computable "60 "(premise_invalid); no proxy window or redefined treatment used"})61 feas = "infeasible"62 reason = (f"The file covers {start} to {end} and contains 0 rows with post_cut = 1, so neither the post-cut "63 "mean change nor the 120-day pre-window (anchored on the first flagged day) can be computed.")64 conclusion = "premise_invalid"65 excl.append(dict(step=3, rule="keep post_cut=1 days and the 120 rows before the first flagged day",66 n_before=len(chg), n_after=0, n_removed=len(chg),67 rationale="no flagged day exists, so the estimation windows are empty"))68else:69 first = int(np.flatnonzero(df["post_cut"].values == 1)[0])70 post = df.loc[df["post_cut"] == 1, "chg_bp"].dropna()71 pre = df["chg_bp"].iloc[max(first - PRE_N, 1):first].dropna()72 if len(pre) < PRE_N:73 issues.append({"issue": f"pre-window has only {len(pre)} changes (<{PRE_N})", "handling": "used as available"})74 excl.append(dict(step=3, rule="keep post_cut=1 days and the 120 rows before the first flagged day",75 n_before=len(chg), n_after=len(pre) + len(post),76 n_removed=len(chg) - len(pre) - len(post), rationale="estimation windows per TASK"))77 if len(pre) >= 2 and len(post) >= 2:78 effect = float(post.mean() - pre.mean())79 v1, v0 = post.var(ddof=1) / len(post), pre.var(ddof=1) / len(pre)80 se = float(np.sqrt(v1 + v0))81 dof = (v1 + v0) ** 2 / (v1 ** 2 / (len(post) - 1) + v0 ** 2 / (len(pre) - 1))82 q = stats.t.ppf(0.975, dof)83 lo, hi = effect - q * se, effect + q * se84 resid = np.concatenate([post - post.mean(), pre - pre.mean()])85 rmse = float(np.sqrt(np.mean(resid ** 2)))86 n_analysed = int(len(pre) + len(post))87 conclusion = "positive_effect" if lo > 0 else "negative_effect" if hi < 0 else "no_detectable_effect"88 else:89 feas, conclusion = "infeasible", "premise_invalid"90 reason = "fewer than two observations in a window"9192results = {93 "task_id": "F6-conf-02", "sample_start": start, "sample_end": end, "n_treated": n_treated,94 "baseline_metric": baseline, "model_metric": rmse, "effect_estimate": effect,95 "effect_ci_low": lo, "effect_ci_high": hi, "conclusion": conclusion,96 "feasibility": feas, "feasibility_reason": reason, "data_issues": issues,97 "n_input_rows": {"gc5y_daily.csv": n_raw}, "n_analysed": n_analysed, "random_seed": None,98}99with open(f"{OUT}/results.json", "w") as f:100 json.dump(results, f, indent=2)101pd.DataFrame(excl).to_csv(f"{OUT}/exclusions.csv", index=False)102103# ---- table 1 ----104df["year"] = df["dt"].dt.year105t1 = df.groupby("year").agg(n_days=("yield_5y", "size"), mean_yield=("yield_5y", "mean"),106 mean_change_bp=("chg_bp", "mean")).reset_index()107t1.round(6).to_csv(f"{OUT}/table1.csv", index=False)108109# ---- figure 1 ----110fig, ax = plt.subplots(figsize=(9, 4))111ax.plot(df["dt"], df["yield_5y"], lw=0.8)112ax.set_title("GoC 5-year benchmark yield (BD.CDN.5YR.DQ.YLD), file sample\n"113 f"{start} to {end}; rows with post_cut=1: {n_treated}")114ax.set_ylabel("Yield (%)"); ax.set_xlabel("Date")115fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120)116117print(json.dumps({k: results[k] for k in ["sample_start", "sample_end", "n_treated", "baseline_metric",118 "effect_estimate", "conclusion", "feasibility"]}, indent=1))119