1"""Replication entry point for F6-conf-01. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3Reads only data/DGS10.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 statsmodels.formula.api as smf9import matplotlib10matplotlib.use("Agg")11import matplotlib.pyplot as plt1213OUT = "out"; os.makedirs(OUT, exist_ok=True)14DAYS = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]15SPLIT = pd.Timestamp("2009-12-31")1617raw = pd.read_csv("data/DGS10.csv", dtype=str, keep_default_na=False)18n_raw = len(raw)19empty = raw["DGS10"].str.strip() == ""20n_missing = int(empty.sum())21d = raw.loc[~empty].copy()22d["date"] = pd.to_datetime(d["observation_date"])23d["y"] = pd.to_numeric(d["DGS10"], errors="raise")24d = d.sort_values("date").reset_index(drop=True)25assert not d["date"].duplicated().any()26d["dy"] = 100.0 * d["y"].diff()27d["gap"] = d["date"].diff().dt.days28ch = d.iloc[1:].copy() # first available observation has no change29ch["wd"] = ch["date"].dt.dayofweek30assert ch["wd"].between(0, 4).all()31n_changes = len(ch)3233excl = pd.DataFrame([34 dict(step=1, rule="drop rows with empty DGS10 (no yield published)", n_before=n_raw,35 n_after=n_raw - n_missing, n_removed=n_missing, rationale="TASK rule: drop empty rows first"),36 dict(step=2, rule="first available observation has no one-day change", n_before=n_raw - n_missing,37 n_after=n_changes, n_removed=n_raw - n_missing - n_changes, rationale="TASK rule: change needs y_{t-1}"),38])39excl.to_csv(f"{OUT}/exclusions.csv", index=False)4041train = ch[ch["date"] <= SPLIT]; test = ch[ch["date"] > SPLIT]42assert train["date"].max() < test["date"].min()4344# (i) classical one-way ANOVA, all changes45groups = [ch.loc[ch.wd == k, "dy"].to_numpy() for k in range(5)]46F, p = stats.f_oneway(*groups)47# robustness (not a required key): Welch-type check via Kruskal-Wallis and Levene, reported in stdout only48kw = stats.kruskal(*groups); lev = stats.levene(*groups)4950# (ii) Monday contrast, Welch two-sample t interval51mon = ch.loc[ch.wd == 0, "dy"].to_numpy(); oth = ch.loc[ch.wd != 0, "dy"].to_numpy()52eff = mon.mean() - oth.mean()53v1, v2 = mon.var(ddof=1) / len(mon), oth.var(ddof=1) / len(oth)54se = np.sqrt(v1 + v2)55df = (v1 + v2) ** 2 / (v1 ** 2 / (len(mon) - 1) + v2 ** 2 / (len(oth) - 1))56tq = stats.t.ppf(0.975, df)57lo, hi = eff - tq * se, eff + tq * se5859# robustness (extra keys): Monday-indicator regression with HC1 and Newey-West(5) SEs; ANOVA excluding 2008-09 to 2009-06 crisis window60ch["mon"] = (ch["wd"] == 0).astype(int)61hc1 = smf.ols("dy ~ mon", ch).fit(cov_type="HC1").conf_int().loc["mon"].tolist()62nw5 = smf.ols("dy ~ mon", ch).fit(cov_type="HAC", cov_kwds={"maxlags": 5}).conf_int().loc["mon"].tolist()63xc = ch[(ch["date"] < "2008-09-01") | (ch["date"] > "2009-06-30")]64p_xc = stats.f_oneway(*[xc.loc[xc.wd == k, "dy"] for k in range(5)]).pvalue6566# (iii) out-of-sample RMSE67wmeans = train.groupby("wd")["dy"].mean()68fc = test["wd"].map(wmeans).to_numpy()69rmse0 = float(np.sqrt(np.mean(test["dy"].to_numpy() ** 2)))70rmse1 = float(np.sqrt(np.mean((test["dy"].to_numpy() - fc) ** 2)))7172if p < 0.05:73 concl = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")74else:75 concl = "no_detectable_effect"7677# table and figure78rows = []79for k, g in enumerate(groups):80 rows.append(dict(weekday=DAYS[k], n=len(g), mean_bp=g.mean(), se_bp=g.std(ddof=1) / np.sqrt(len(g))))81tab = pd.DataFrame(rows)82tab.to_csv(f"{OUT}/table1.csv", index=False)83fig, ax = plt.subplots(figsize=(6, 4))84ci = [stats.t.ppf(0.975, r["n"] - 1) * r["se_bp"] for r in rows]85ax.errorbar(range(5), tab["mean_bp"], yerr=ci, fmt="o", capsize=5)86ax.axhline(0, color="grey", lw=0.8)87ax.set_xticks(range(5)); ax.set_xticklabels([s[:3] for s in DAYS])88ax.set_ylabel("Mean one-day change in DGS10 (bp)")89ax.set_title("10y UST yield: mean daily change by weekday, 2000-2019\n(95%% CI, all %d changes)" % n_changes)90fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120); plt.close(fig)9192gap_counts = ch["gap"].value_counts().sort_index()93res = {94 "task_id": "F6-conf-01",95 "feasibility": "complete",96 "feasibility_reason": None,97 "data_issues": [98 {"issue": f"{n_missing} rows with empty DGS10 (US market holidays incl. 2001-09-11/12)",99 "handling": "dropped before differencing, per TASK rules"},100 {"issue": f"{int((ch.gap > 3).sum() + ((ch.gap == 2)).sum())} changes span a holiday (calendar gap 2 or 4 days), i.e. cover more than one trading session",101 "handling": "kept and attributed to weekday of t as the TASK prescribes; noted as a limitation"},102 {"issue": "extreme changes (e.g. -51 bp on 2009-03-18, -27 bp on 2008-09-15) during the 2008-09 crisis",103 "handling": "kept (genuine market events); no outlier trimming; discussed as limitation"},104 {"issue": "yields quoted to 0.01 percentage point, so changes are discretised at 1 bp",105 "handling": "none needed; noted"},106 ],107 "n_input_rows": {"DGS10.csv": int(n_raw)},108 "n_analysed": int(n_changes),109 "random_seed": None,110 "n_missing_rows": n_missing,111 "n_train": int(len(train)),112 "n_test": int(len(test)),113 "anova_f_stat": float(F),114 "anova_p_value": float(p),115 "effect_estimate": float(eff),116 "effect_ci_low": float(lo),117 "effect_ci_high": float(hi),118 "baseline_metric": rmse0,119 "model_metric": rmse1,120 "conclusion": concl,121 "effect_se_welch": float(se),122 "effect_df_welch": float(df),123 "train_weekday_means_bp": {DAYS[k]: float(v) for k, v in wmeans.items()},124 "rmse_ratio_model_over_baseline": rmse1 / rmse0,125 "robust_kruskal_wallis_p": float(kw.pvalue),126 "robust_levene_p": float(lev.pvalue),127 "robust_effect_ci_hc1": [float(v) for v in hc1],128 "robust_effect_ci_nw5": [float(v) for v in nw5],129 "robust_anova_p_excl_crisis_2008_09_to_2009_06": float(p_xc),130 "n_excl_crisis": int(len(xc)),131}132with open(f"{OUT}/results.json", "w") as f:133 json.dump(res, f, indent=2)134print(json.dumps({k: res[k] for k in ["n_train", "n_test", "anova_f_stat", "anova_p_value", "effect_estimate",135 "effect_ci_low", "effect_ci_high", "baseline_metric", "model_metric", "conclusion"]}, indent=1))136