1"""Replication entry point for F6-conf-01. Run from workspace root."""2import json, os3import numpy as np, pandas as pd4from scipy import stats5import matplotlib6matplotlib.use("Agg")7import matplotlib.pyplot as plt89OUT = "out"10os.makedirs(OUT, exist_ok=True)11raw = pd.read_csv("data/DGS10.csv", dtype=str, keep_default_na=False)12n_raw = len(raw)13missing = raw["DGS10"].str.strip() == ""14n_missing = int(missing.sum())15d = raw[~missing].copy()16d["date"] = pd.to_datetime(d["observation_date"])17d["y"] = d["DGS10"].astype(float)18d = d.sort_values("date").reset_index(drop=True)19gapdays = d["date"].diff().dt.days20d["dy"] = 100 * d["y"].diff()21ch = d.iloc[1:].copy()22ch["dy"] = ch["dy"].round(10) # remove float noise from 2-decimal inputs23ch["wd"] = ch["date"].dt.dayofweek24names = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]25n_changes = len(ch)2627# exclusion ledger28ex = pd.DataFrame([29 ["1", "drop rows with empty DGS10", n_raw, n_raw - n_missing, n_missing, "no yield published (market holiday)"],30 ["2", "first available observation has no change", n_raw - n_missing, n_changes, 1, "change undefined"],31], columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"])32ex.to_csv(f"{OUT}/exclusions.csv", index=False)3334train = ch[ch["date"] <= "2009-12-31"]35test = ch[ch["date"] >= "2010-01-01"]3637# ANOVA (classical, equal variance), all changes38groups = [ch.loc[ch.wd == k, "dy"].values for k in range(5)]39F, p = stats.f_oneway(*groups)4041# Monday contrast, Welch42mon = ch.loc[ch.wd == 0, "dy"].values43oth = ch.loc[ch.wd != 0, "dy"].values44est = mon.mean() - oth.mean()45v1, v2 = mon.var(ddof=1) / len(mon), oth.var(ddof=1) / len(oth)46se = np.sqrt(v1 + v2)47df = (v1 + v2) ** 2 / (v1 ** 2 / (len(mon) - 1) + v2 ** 2 / (len(oth) - 1))48tq = stats.t.ppf(0.975, df)49lo, hi = est - tq * se, est + tq * se5051# Out-of-sample52wd_means = train.groupby("wd")["dy"].mean()53pred = test["wd"].map(wd_means).values54rmse_base = float(np.sqrt(np.mean(test["dy"].values ** 2)))55rmse_model = float(np.sqrt(np.mean((test["dy"].values - pred) ** 2)))5657if p >= 0.05:58 concl = "no_detectable_effect"59elif lo > 0:60 concl = "positive_effect"61elif hi < 0:62 concl = "negative_effect"63else:64 concl = "no_detectable_effect"6566# Table and figure67tab = pd.DataFrame({68 "weekday": names,69 "n": [len(g) for g in groups],70 "mean_bp": [g.mean() for g in groups],71 "se_bp": [g.std(ddof=1) / np.sqrt(len(g)) for g in groups],72})73tab.to_csv(f"{OUT}/table1.csv", index=False)74fig, ax = plt.subplots(figsize=(6, 4))75crit = [stats.t.ppf(0.975, n - 1) for n in tab["n"]]76ax.errorbar(range(5), tab["mean_bp"], yerr=np.array(crit) * tab["se_bp"], fmt="o", capsize=4)77ax.axhline(0, color="grey", lw=0.8)78ax.set_xticks(range(5)); ax.set_xticklabels([n[:3] for n in names])79ax.set_ylabel("Mean daily change in 10y yield (bp)")80ax.set_title("DGS10 daily change by weekday, 2000-2019 (95% CI)")81fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120); plt.close(fig)8283n_multiday = int(((gapdays.iloc[1:] == 2) | ((gapdays.iloc[1:] >= 3) & (d["date"].iloc[1:].dt.dayofweek != 0))84 | (gapdays.iloc[1:] >= 4)).sum())85res = {86 "task_id": "F6-conf-01",87 "feasibility": "complete",88 "feasibility_reason": None,89 "data_issues": [90 {"issue": f"{n_missing} rows with empty DGS10 (US market holidays, incl. 2001-09-11/12)",91 "handling": "dropped before differencing, as required by TASK.md"},92 {"issue": f"{n_multiday} changes span more than one trading day because the previous weekday was a holiday",93 "handling": "kept and attributed to weekday of t, per TASK.md rules; noted as limitation"},94 ],95 "n_input_rows": {"DGS10.csv": int(n_raw)},96 "n_analysed": int(n_changes),97 "random_seed": None,98 "n_missing_rows": n_missing,99 "n_train": int(len(train)),100 "n_test": int(len(test)),101 "anova_f_stat": float(F),102 "anova_p_value": float(p),103 "effect_estimate": float(est),104 "effect_ci_low": float(lo),105 "effect_ci_high": float(hi),106 "baseline_metric": rmse_base,107 "model_metric": rmse_model,108 "conclusion": concl,109 "effect_se_welch": float(se),110 "effect_df_welch": float(df),111 "train_weekday_means_bp": {names[k]: float(v) for k, v in wd_means.items()},112 "rmse_ratio_model_over_baseline": rmse_model / rmse_base,113 "test_mean_bp": float(test["dy"].mean()),114}115with open(f"{OUT}/results.json", "w") as f:116 json.dump(res, f, indent=2)117print(json.dumps({k: res[k] for k in ["n_train", "n_test", "anova_f_stat", "anova_p_value", "effect_estimate",118 "effect_ci_low", "effect_ci_high", "baseline_metric", "model_metric", "conclusion"]}, indent=1))119