1"""F6-conf-01: weekday effects in daily changes of the 10-year 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 matplotlib9matplotlib.use("Agg")10import matplotlib.pyplot as plt1112OUT = "out"; os.makedirs(OUT, exist_ok=True)13DAYS = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]1415raw = pd.read_csv("data/DGS10.csv", dtype=str, keep_default_na=False)16n_raw = len(raw)17raw["date"] = pd.to_datetime(raw["observation_date"], format="%Y-%m-%d")18assert raw["date"].is_monotonic_increasing and not raw["date"].duplicated().any()19assert (raw["date"].dt.dayofweek < 5).all()20empty = raw["DGS10"].str.strip() == ""21n_missing = int(empty.sum())22obs = raw.loc[~empty].copy()23obs["y"] = obs["DGS10"].astype(float)24obs = obs.reset_index(drop=True)2526# One-day change between consecutive available observations, dated t27obs["dy"] = 100.0 * obs["y"].diff()28obs["gap_days"] = obs["date"].diff().dt.days29ch = obs.iloc[1:].copy() # first available observation has no change30ch["dy"] = ch["dy"].round(10) # remove float noise from 2-decimal percent data31ch["wd"] = ch["date"].dt.dayofweek32n_changes = len(ch)3334excl = pd.DataFrame([35 ["1", "drop rows with empty DGS10 (market holidays)", n_raw, n_raw - n_missing, n_missing,36 "TASK rule: no yield published; dropped before differencing"],37 ["2", "first available observation has no change", n_raw - n_missing, n_changes, 1,38 "Δy undefined for first observation"],39], columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"])40excl.to_csv(f"{OUT}/exclusions.csv", index=False)4142train = ch[ch["date"] <= "2009-12-31"]; test = ch[ch["date"] >= "2010-01-01"]43assert len(train) + len(test) == n_changes4445# (i) one-way ANOVA (equal variance), all changes46groups = [ch.loc[ch.wd == k, "dy"].to_numpy() for k in range(5)]47F, p = stats.f_oneway(*groups)4849# Table 150tab = pd.DataFrame({"weekday": DAYS, "n": [len(g) for g in groups],51 "mean_bp": [g.mean() for g in groups],52 "se_bp": [g.std(ddof=1) / np.sqrt(len(g)) for g in groups]})53tab.to_csv(f"{OUT}/table1.csv", index=False, float_format="%.6f")5455# (ii) Monday contrast with Welch 95% CI56mon, oth = groups[0], np.concatenate(groups[1:])57est = mon.mean() - oth.mean()58v1, v2 = mon.var(ddof=1) / len(mon), oth.var(ddof=1) / len(oth)59se = np.sqrt(v1 + v2)60df_w = (v1 + v2) ** 2 / (v1 ** 2 / (len(mon) - 1) + v2 ** 2 / (len(oth) - 1))61tq = stats.t.ppf(0.975, df_w)62lo, hi = est - tq * se, est + tq * se63welch_p = 2 * stats.t.sf(abs(est / se), df_w)6465# (iii) out-of-sample RMSE66wmeans = train.groupby("wd")["dy"].mean()67pred = test["wd"].map(wmeans).to_numpy()68rmse0 = float(np.sqrt(np.mean(test["dy"].to_numpy() ** 2)))69rmse1 = float(np.sqrt(np.mean((test["dy"].to_numpy() - pred) ** 2)))70# Diebold-Mariano-type check on squared-error loss differential (HAC not needed for h=1; iid SE)71d = (test["dy"].to_numpy() - pred) ** 2 - test["dy"].to_numpy() ** 272dm_t = d.mean() / (d.std(ddof=1) / np.sqrt(len(d)))73dm_p = 2 * stats.norm.sf(abs(dm_t))7475if p < 0.05:76 conclusion = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")77else:78 conclusion = "no_detectable_effect"7980gap_long = ch[ch["gap_days"] > 3]81res = {82 "task_id": "F6-conf-01",83 "feasibility": "complete",84 "feasibility_reason": None,85 "data_issues": [86 {"issue": f"{n_missing} rows with empty DGS10 field (US market holidays)",87 "handling": "dropped before differencing, per TASK rule"},88 {"issue": f"{len(gap_long)} changes span more than 3 calendar days (holiday gaps, e.g. after Monday holidays or 2001-09-11 closure)",89 "handling": "kept as one-day changes between consecutive available observations, attributed to weekday of t, per TASK rule"},90 {"issue": "heavy-tailed changes (max |dy| = %.0f bp on %s)" % (ch.dy.abs().max(), ch.loc[ch.dy.abs().idxmax(), "observation_date"]),91 "handling": "kept: large moves are plausible market events (the largest coincides with the 2009-03-18 FOMC announcement); no outlier removal; Welch CI used for the contrast"},92 ],93 "n_input_rows": {"DGS10.csv": int(n_raw)},94 "n_analysed": int(n_changes),95 "random_seed": None,96 "n_missing_rows": n_missing,97 "n_train": int(len(train)),98 "n_test": int(len(test)),99 "anova_f_stat": float(F),100 "anova_p_value": float(p),101 "effect_estimate": float(est),102 "effect_ci_low": float(lo),103 "effect_ci_high": float(hi),104 "baseline_metric": rmse0,105 "model_metric": rmse1,106 "conclusion": conclusion,107 # supplementary108 "effect_se": float(se),109 "effect_welch_df": float(df_w),110 "effect_welch_p_value": float(welch_p),111 "oos_mse_diff_mean": float(d.mean()),112 "oos_dm_t_stat": float(dm_t),113 "oos_dm_p_value": float(dm_p),114 "train_weekday_means_bp": {DAYS[k]: float(wmeans[k]) for k in range(5)},115 "rmse_ratio_model_to_baseline": rmse1 / rmse0,116}117with open(f"{OUT}/results.json", "w") as f:118 json.dump(res, f, indent=2)119120# Figure 1121ci = stats.t.ppf(0.975, tab["n"] - 1) * tab["se_bp"]122fig, ax = plt.subplots(figsize=(6, 4))123ax.errorbar(range(5), tab["mean_bp"], yerr=ci, fmt="o", capsize=5)124ax.axhline(0, color="grey", lw=0.8)125ax.set_xticks(range(5)); ax.set_xticklabels(DAYS)126ax.set_ylabel("Mean one-day change in DGS10 (bp)")127ax.set_title(f"10y Treasury yield change by weekday, 2000-2019 (n={n_changes})\nANOVA F={F:.2f}, p={p:.3f}; bars = 95% CI")128fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120)129print(json.dumps({k: res[k] for k in ["n_analysed", "anova_f_stat", "anova_p_value", "effect_estimate",130 "effect_ci_low", "effect_ci_high", "baseline_metric", "model_metric", "conclusion"]}, indent=1))131