1"""F6-conf-01: weekday effects in daily changes of the 10y Treasury yield (DGS10).2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/DGS10.csv; writes out/results.json, out/table1.csv, out/figure1.png, out/exclusions.csv.4No randomness is used (random_seed = null)."""5import json, os6import numpy as np, pandas as pd7from scipy import stats8import statsmodels.api as sm9import matplotlib10matplotlib.use("Agg")11import matplotlib.pyplot as plt1213OUT = "out"; os.makedirs(OUT, exist_ok=True)14DAYS = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]1516raw = pd.read_csv("data/DGS10.csv", dtype=str, keep_default_na=False)17n_raw = len(raw)18empty = raw["DGS10"].str.strip() == ""19n_missing = int(empty.sum())20d = raw[~empty].copy()21d["date"] = pd.to_datetime(d["observation_date"], format="%Y-%m-%d")22d["y"] = pd.to_numeric(d["DGS10"], errors="raise")23assert d["date"].is_monotonic_increasing and not d["date"].duplicated().any()24d["dy"] = 100.0 * d["y"].diff()25d["dy"] = d["dy"].round(8) # remove float noise; data have 1 bp resolution26ch = d.dropna(subset=["dy"]).copy()27ch["wd"] = ch["date"].dt.dayofweek28assert ch["wd"].between(0, 4).all()29ch["weekday"] = ch["wd"].map(dict(enumerate(DAYS)))3031pd.DataFrame([32 dict(step=1, rule="drop rows with empty DGS10 field (holidays)", n_before=n_raw,33 n_after=n_raw - n_missing, n_removed=n_missing, rationale="TASK rule: no published yield"),34 dict(step=2, rule="first available observation has no change", n_before=n_raw - n_missing,35 n_after=len(ch), n_removed=n_raw - n_missing - len(ch), rationale="TASK rule: change needs y_{t-1}"),36]).to_csv(f"{OUT}/exclusions.csv", index=False)3738train = ch[ch["date"] <= "2009-12-31"]; test = ch[ch["date"] >= "2010-01-01"]39assert len(train) + len(test) == len(ch)4041# (i) classical one-way ANOVA, all changes42groups = [ch.loc[ch["wd"] == k, "dy"].values for k in range(5)]43F, p = stats.f_oneway(*groups)44kw_H, kw_p = stats.kruskal(*groups) # sensitivity (rank-based)45lev_W, lev_p = stats.levene(*groups, center="median") # equal-variance check46from statsmodels.stats.oneway import anova_oneway47wa = anova_oneway(groups, use_var="unequal") # Welch ANOVA sensitivity4849# (ii) Monday contrast, Welch two-sample interval (primary)50mon = ch.loc[ch["wd"] == 0, "dy"].values; oth = ch.loc[ch["wd"] != 0, "dy"].values51est = mon.mean() - oth.mean()52v1, v2 = mon.var(ddof=1) / len(mon), oth.var(ddof=1) / len(oth)53se = np.sqrt(v1 + v2)54df_w = (v1 + v2) ** 2 / (v1 ** 2 / (len(mon) - 1) + v2 ** 2 / (len(oth) - 1))55tq = stats.t.ppf(0.975, df_w)56lo, hi = est - tq * se, est + tq * se57welch_p = 2 * stats.t.sf(abs(est / se), df_w)58# sensitivity: OLS on Monday dummy with HAC (Newey-West, 5 lags) SE59X = sm.add_constant((ch["wd"] == 0).astype(float).values)60hac = sm.OLS(ch["dy"].values, X).fit(cov_type="HAC", cov_kwds={"maxlags": 5})6162# (iii) out-of-sample63wmeans = train.groupby("wd")["dy"].mean()64f_model = test["wd"].map(wmeans).values65e0 = test["dy"].values; e1 = e0 - f_model66rmse0 = float(np.sqrt(np.mean(e0 ** 2))); rmse1 = float(np.sqrt(np.mean(e1 ** 2)))67e_unc = e0 - train["dy"].mean()68rmse_unc = float(np.sqrt(np.mean(e_unc ** 2)))69# Diebold-Mariano on squared-error loss differential (model - baseline), HAC 5 lags70dl = e1 ** 2 - e0 ** 271dm = sm.OLS(dl, np.ones_like(dl)).fit(cov_type="HAC", cov_kwds={"maxlags": 5})7273if p >= 0.05:74 concl = "no_detectable_effect"75elif lo > 0:76 concl = "positive_effect"77elif hi < 0:78 concl = "negative_effect"79else:80 concl = "no_detectable_effect"8182# table + figure83tab = ch.groupby("wd")["dy"].agg(n="count", mean_bp="mean", sd="std").reindex(range(5))84tab["se_bp"] = tab["sd"] / np.sqrt(tab["n"])85tab.insert(0, "weekday", DAYS)86tab[["weekday", "n", "mean_bp", "se_bp"]].to_csv(f"{OUT}/table1.csv", index=False, float_format="%.6f")87crit = stats.t.ppf(0.975, tab["n"] - 1)88fig, ax = plt.subplots(figsize=(6, 4))89ax.errorbar(range(5), tab["mean_bp"], yerr=crit * tab["se_bp"], fmt="o", capsize=5, color="C0")90ax.axhline(0, color="grey", lw=0.8, ls="--")91ax.set_xticks(range(5)); ax.set_xticklabels([s[:3] for s in DAYS])92ax.set_ylabel("Mean one-day change in DGS10 (bp)")93ax.set_title(f"10y Treasury yield changes by weekday, 2000-2019 (n={len(ch)})\n"94 f"95% CI; one-way ANOVA F={F:.2f}, p={p:.3f}", fontsize=9)95fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=150); plt.close(fig)9697r = lambda x: float(round(float(x), 10))98res = {99 "task_id": "F6-conf-01", "feasibility": "complete", "feasibility_reason": None,100 "data_issues": [101 {"issue": f"{n_missing} rows with empty DGS10 field (US market holidays)",102 "handling": "dropped before differencing, per TASK rule"},103 {"issue": "changes following a holiday span more than one trading day (e.g. Fri->Tue after Monday holidays)",104 "handling": "kept and attributed to weekday of t, per TASK rule; noted as limitation"},105 {"issue": "extreme changes (min -51 bp on 2009-03-18, QE1 announcement; max +25 bp); heavy tails",106 "handling": "kept (genuine market events); rank-based Kruskal-Wallis reported as sensitivity"},107 {"issue": "yields recorded to 0.01 pct, so changes are discretised at 1 bp",108 "handling": "none needed; noted"},109 ],110 "n_input_rows": {"DGS10.csv": int(n_raw)},111 "n_analysed": int(len(ch)), "random_seed": None,112 "n_missing_rows": n_missing, "n_train": int(len(train)), "n_test": int(len(test)),113 "anova_f_stat": r(F), "anova_p_value": r(p),114 "effect_estimate": r(est), "effect_ci_low": r(lo), "effect_ci_high": r(hi),115 "baseline_metric": r(rmse0), "model_metric": r(rmse1),116 "conclusion": concl,117 # supplementary118 "effect_se_welch": r(se), "effect_df_welch": r(df_w), "effect_p_welch": r(welch_p),119 "effect_ci_method": "Welch two-sample t interval (Monday vs Tue-Fri), all changes",120 "sens_hac_effect_se": r(hac.bse[1]), "sens_hac_effect_p": r(hac.pvalues[1]),121 "sens_kruskal_H": r(kw_H), "sens_kruskal_p": r(kw_p),122 "sens_welch_anova_F": r(wa.statistic), "sens_welch_anova_p": r(wa.pvalue),123 "sens_levene_W": r(lev_W), "sens_levene_p": r(lev_p),124 "oos_rmse_ratio_model_over_baseline": r(rmse1 / rmse0),125 "oos_rmse_unconditional_train_mean": r(rmse_unc),126 "oos_dm_mean_loss_diff": r(dm.params[0]), "oos_dm_t": r(dm.tvalues[0]), "oos_dm_p": r(dm.pvalues[0]),127 "train_weekday_means_bp": {DAYS[k]: r(v) for k, v in wmeans.items()},128 "n_changes_total": int(len(ch)),129}130with open(f"{OUT}/results.json", "w") as f:131 json.dump(res, f, indent=2)132print(json.dumps({k: res[k] for k in ["n_analysed", "anova_f_stat", "anova_p_value", "effect_estimate",133 "effect_ci_low", "effect_ci_high", "baseline_metric", "model_metric", "conclusion"]}, indent=1))134