1"""Replication entry point for F6-conf-01. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py"""3import json, os4import numpy as np, pandas as pd5from scipy import stats6import matplotlib7matplotlib.use("Agg")8import matplotlib.pyplot as plt910SEED = None # no randomness used11os.makedirs("out", exist_ok=True)12raw = pd.read_csv("data/DGS10.csv", dtype=str, keep_default_na=False)13n_input = len(raw)14assert not raw.observation_date.duplicated().any()15n_missing = int((raw.DGS10.str.strip() == "").sum())16d = raw[raw.DGS10.str.strip() != ""].copy()17d["date"] = pd.to_datetime(d.observation_date)18d["y"] = d.DGS10.astype(float)19d = d.sort_values("date").reset_index(drop=True)20d["dy"] = 100 * d.y.diff()21d["gap_days"] = d.date.diff().dt.days22d = d.dropna(subset=["dy"]).copy()23d["dy"] = d.dy.round(10)24d["wd"] = d.date.dt.dayofweek25names = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]26assert d.wd.isin(range(5)).all()2728train = d[d.date <= "2009-12-31"]; test = d[d.date >= "2010-01-01"]2930# (i) one-way ANOVA, all changes31groups = [d.dy[d.wd == k].values for k in range(5)]32F, p = stats.f_oneway(*groups)3334# (ii) Monday contrast, Welch two-sample interval35mon = d.dy[d.wd == 0].values; oth = d.dy[d.wd != 0].values36eff = mon.mean() - oth.mean()37v1, v2 = mon.var(ddof=1) / len(mon), oth.var(ddof=1) / len(oth)38se = np.sqrt(v1 + v2)39df = (v1 + v2) ** 2 / (v1 ** 2 / (len(mon) - 1) + v2 ** 2 / (len(oth) - 1))40tq = stats.t.ppf(0.975, df)41lo, hi = eff - tq * se, eff + tq * se4243# (iii) out-of-sample44wd_means = train.groupby("wd").dy.mean()45pred = test.wd.map(wd_means).values46rmse0 = float(np.sqrt(np.mean(test.dy.values ** 2)))47rmse1 = float(np.sqrt(np.mean((test.dy.values - pred) ** 2)))48# Diebold-Mariano (squared loss, h=1) as supplementary49dl = (test.dy.values - pred) ** 2 - test.dy.values ** 250dm = dl.mean() / (dl.std(ddof=1) / np.sqrt(len(dl)))51dm_p = 2 * stats.norm.sf(abs(dm))52# supplementary: train-only mean (unconditional) forecast53rmse_mean = float(np.sqrt(np.mean((test.dy.values - train.dy.mean()) ** 2)))5455if p < 0.05:56 concl = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")57else:58 concl = "no_detectable_effect"5960# table + figure61tab = pd.DataFrame({"weekday": names,62 "n": [len(g) for g in groups],63 "mean_bp": [g.mean() for g in groups],64 "se_bp": [g.std(ddof=1) / np.sqrt(len(g)) for g in groups]})65tab.to_csv("out/table1.csv", index=False, float_format="%.6f")66fig, ax = plt.subplots(figsize=(6, 4))67ci = [stats.t.ppf(0.975, n - 1) * s for n, s in zip(tab.n, tab.se_bp)]68ax.errorbar(range(5), tab.mean_bp, yerr=ci, fmt="o", capsize=5)69ax.axhline(0, color="grey", lw=0.8)70ax.set_xticks(range(5)); ax.set_xticklabels(names)71ax.set_ylabel("Mean one-day change in 10y yield (bp)")72ax.set_title(f"DGS10 daily changes by weekday, 2000-2019 (n={len(d)})\n95% CIs; ANOVA F={F:.2f}, p={p:.3f}")73fig.tight_layout(); fig.savefig("out/figure1.png", dpi=150)7475post_hol = int(((d.gap_days > 1) & ~((d.wd == 0) & (d.gap_days == 3))).sum())76res = {77 "task_id": "F6-conf-01",78 "n_missing_rows": n_missing,79 "n_train": int(len(train)), "n_test": int(len(test)),80 "anova_f_stat": float(F), "anova_p_value": float(p),81 "effect_estimate": float(eff), "effect_ci_low": float(lo), "effect_ci_high": float(hi),82 "baseline_metric": rmse0, "model_metric": rmse1,83 "conclusion": concl,84 "feasibility": "complete", "feasibility_reason": None,85 "data_issues": [86 {"issue": f"{n_missing} rows with empty DGS10 field (US market holidays / no publication)",87 "handling": "Dropped before differencing, as prescribed; changes computed between consecutive available observations."},88 {"issue": f"{post_hol} changes (excluding ordinary Fri->Mon weekends) span a holiday gap and thus cover more than one business day",89 "handling": "Kept and attributed to the weekday of t, as prescribed by TASK.md; noted as a limitation (inflates variance, mostly for Tuesdays after Monday holidays)."},90 ],91 "n_input_rows": {"DGS10.csv": n_input},92 "n_analysed": int(len(d)),93 "random_seed": SEED,94 "supplementary": {95 "monday_contrast_method": "Welch two-sample t interval (unequal variances)",96 "welch_se_bp": float(se), "welch_df": float(df),97 "welch_p_value": float(2 * stats.t.sf(abs(eff / se), df)),98 "train_weekday_means_bp": {names[k]: float(v) for k, v in wd_means.items()},99 "rmse_train_unconditional_mean_bp": rmse_mean,100 "dm_stat_model_minus_zero_sqloss": float(dm), "dm_p_value": float(dm_p),101 "levene_p_value": float(stats.levene(*groups).pvalue),102 "kruskal_p_value": float(stats.kruskal(*groups).pvalue),103 },104}105with open("out/results.json", "w") as f:106 json.dump(res, f, indent=2)107print(json.dumps(res, indent=2)); print(tab)108