1"""Day-of-week effect in daily 10y Treasury yield changes (F6-conf-01).2Run from workspace root: ~/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 used11OUT = "out"12raw = pd.read_csv("data/DGS10.csv", dtype={"DGS10": str}, keep_default_na=False)13n_input = len(raw)14issues = []15raw["observation_date"] = pd.to_datetime(raw["observation_date"])16val = raw["DGS10"].str.strip()17missing = (val == "") | (val == ".")18n_missing = int(missing.sum())19issues.append({"issue": f"{n_missing} rows with empty DGS10 field (US market holidays / no publication)",20 "handling": "Dropped before computing changes, per TASK rules; the change after a holiday spans the gap and is attributed to the weekday of its end date."})21df = raw.loc[~missing].copy()22df["y"] = pd.to_numeric(df["DGS10"], errors="coerce")23bad = int(df["y"].isna().sum())24if bad:25 issues.append({"issue": f"{bad} non-numeric DGS10 values", "handling": "Dropped."})26df = df.dropna(subset=["y"])27ndup = int(df["observation_date"].duplicated().sum())28nwkend = int((df["observation_date"].dt.dayofweek >= 5).sum())29if ndup or nwkend:30 issues.append({"issue": f"{ndup} duplicate dates, {nwkend} weekend dates", "handling": "Reported only; none expected."})31df = df.sort_values("observation_date").reset_index(drop=True)32df["dy"] = 100 * df["y"].diff()33gap = df["observation_date"].diff().dt.days34n_multi = int(((gap > 1) & (df["observation_date"].dt.dayofweek != 0)).sum() + ((gap > 3) & (df["observation_date"].dt.dayofweek == 0)).sum())35issues.append({"issue": f"{n_multi} changes span a holiday (more calendar/business days than usual), e.g. Tuesday after a Monday holiday",36 "handling": "Kept as defined by the TASK (consecutive available observations); noted as a limitation."})37issues.append({"issue": "Yields are quoted to 2 decimals, so changes are discrete in 1 bp steps",38 "handling": "No action; immaterial for means/ANOVA."})39ch = df.dropna(subset=["dy"]).copy()40ch["wd"] = ch["observation_date"].dt.dayofweek41names = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]4243# ANOVA (all changes)44groups = [ch.loc[ch.wd == d, "dy"].values for d in range(5)]45F, p = stats.f_oneway(*groups)4647# Monday contrast: Welch two-sample48mon, oth = groups[0], ch.loc[ch.wd != 0, "dy"].values49eff = mon.mean() - oth.mean()50v1, v2 = mon.var(ddof=1) / len(mon), oth.var(ddof=1) / len(oth)51se = np.sqrt(v1 + v2)52dfw = (v1 + v2) ** 2 / (v1**2 / (len(mon) - 1) + v2**2 / (len(oth) - 1))53t = stats.t.ppf(0.975, dfw)54lo, hi = eff - t * se, eff + t * se55welch_p = 2 * stats.t.sf(abs(eff / se), dfw)5657# Out-of-sample58train = ch[ch.observation_date <= "2009-12-31"]59test = ch[ch.observation_date >= "2010-01-01"]60wmeans = train.groupby("wd")["dy"].mean()61pred = test["wd"].map(wmeans).values62rmse_base = float(np.sqrt(np.mean(test["dy"].values ** 2)))63rmse_model = float(np.sqrt(np.mean((test["dy"].values - pred) ** 2)))64# Diebold-Mariano (squared loss, simple HAC-free since 1-step) for descriptive purposes65d = test["dy"].values ** 2 - (test["dy"].values - pred) ** 266dm_t = d.mean() / (d.std(ddof=1) / np.sqrt(len(d)))67dm_p = 2 * stats.norm.sf(abs(dm_t))6869if p < 0.05:70 concl = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")71else:72 concl = "no_detectable_effect"7374# Table and figure75tab = pd.DataFrame({"weekday": names, "n": [len(g) for g in groups],76 "mean_bp": [g.mean() for g in groups],77 "se_bp": [g.std(ddof=1) / np.sqrt(len(g)) for g in groups]})78tab.to_csv(os.path.join(OUT, "table1.csv"), index=False)79fig, ax = plt.subplots(figsize=(6, 4))80ci = [stats.t.ppf(0.975, n - 1) * s for n, s in zip(tab.n, tab.se_bp)]81ax.errorbar(range(5), tab.mean_bp, yerr=ci, fmt="o", capsize=5)82ax.axhline(0, color="grey", lw=0.8)83ax.set_xticks(range(5)); ax.set_xticklabels(names)84ax.set_ylabel("Mean daily change in 10y yield (bp)")85ax.set_title("DGS10 daily change by weekday, 2000-2019 (95% CI)")86fig.tight_layout(); fig.savefig(os.path.join(OUT, "figure1.png"), dpi=150)8788res = {89 "task_id": "F6-conf-01",90 "n_missing_rows": n_missing,91 "n_train": int(len(train)), "n_test": int(len(test)),92 "anova_f_stat": float(F), "anova_p_value": float(p),93 "effect_estimate": float(eff), "effect_ci_low": float(lo), "effect_ci_high": float(hi),94 "baseline_metric": rmse_base, "model_metric": rmse_model,95 "conclusion": concl, "feasibility": "complete", "feasibility_reason": None,96 "data_issues": issues, "n_input_rows": {"DGS10.csv": int(n_input)},97 "n_analysed": int(len(ch)), "random_seed": SEED,98 "supplementary": {"monday_contrast_method": "Welch two-sample t interval",99 "welch_se_bp": float(se), "welch_df": float(dfw), "welch_p_value": float(welch_p),100 "dm_stat_squared_loss": float(dm_t), "dm_p_value": float(dm_p),101 "train_weekday_means_bp": {names[k]: float(v) for k, v in wmeans.items()},102 "sd_all_changes_bp": float(ch.dy.std(ddof=1))},103}104with open(os.path.join(OUT, "results.json"), "w") as f:105 json.dump(res, f, indent=2)106print(json.dumps(res, indent=2)); print(tab)107