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.py"""3import json, hashlib, os4import numpy as np, pandas as pd5from scipy import stats6import matplotlib7matplotlib.use("Agg")8import matplotlib.pyplot as plt910SEED = None # no randomness used11DATA = "data/DGS10.csv"; OUT = "out"12os.makedirs(OUT, exist_ok=True)13DAYS = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]1415raw = pd.read_csv(DATA, dtype=str, keep_default_na=False)16n_raw = len(raw)17sha = hashlib.sha256(open(DATA, "rb").read()).hexdigest()18raw["date"] = pd.to_datetime(raw["observation_date"], format="%Y-%m-%d")19empty = raw["DGS10"].str.strip() == ""20n_missing = int(empty.sum())21df = raw.loc[~empty].copy()22df["y"] = pd.to_numeric(df["DGS10"], errors="raise")23df = df.sort_values("date").reset_index(drop=True)2425# audit facts26data_issues = []27n_dup = int(raw["observation_date"].duplicated().sum())28n_weekend = int((raw["date"].dt.dayofweek >= 5).sum())29data_issues.append({"issue": f"{n_missing} rows with empty DGS10 field (US market holidays)",30 "handling": "dropped before differencing, per task rule; change after a holiday spans the gap"})31if n_dup: data_issues.append({"issue": f"{n_dup} duplicated dates", "handling": "reported only"})32if n_weekend: data_issues.append({"issue": f"{n_weekend} weekend dates", "handling": "reported only"})3334df["dy"] = 100.0 * df["y"].diff()35df["dy"] = df["dy"].round(10) # remove float noise from percent->bp conversion36ch = df.iloc[1:].copy() # first available obs has no change37ch["wd"] = ch["date"].dt.day_name()38assert set(ch["wd"]) <= set(DAYS)39train = ch[ch["date"] <= "2009-12-31"]; test = ch[ch["date"] >= "2010-01-01"]4041excl = pd.DataFrame([42 ["1", "drop rows with empty DGS10", n_raw, n_raw - n_missing, n_missing, "no published yield (holiday); task rule"],43 ["2", "first available observation has no change", n_raw - n_missing, len(ch), 1, "differencing"],44], columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"])45excl.to_csv(f"{OUT}/exclusions.csv", index=False)4647# (i) ANOVA, all changes48groups = [ch.loc[ch.wd == d, "dy"].values for d in DAYS]49F, p = stats.f_oneway(*groups)5051# table52tab = pd.DataFrame({"weekday": DAYS, "n": [len(g) for g in groups],53 "mean_bp": [g.mean() for g in groups],54 "se_bp": [g.std(ddof=1) / np.sqrt(len(g)) for g in groups]})55tab.to_csv(f"{OUT}/table1.csv", index=False)5657# (ii) Monday contrast, Welch two-sample58mon = ch.loc[ch.wd == "Monday", "dy"].values; oth = ch.loc[ch.wd != "Monday", "dy"].values59eff = mon.mean() - oth.mean()60v1, v2 = mon.var(ddof=1) / len(mon), oth.var(ddof=1) / len(oth)61se = np.sqrt(v1 + v2)62dfw = (v1 + v2) ** 2 / (v1 ** 2 / (len(mon) - 1) + v2 ** 2 / (len(oth) - 1))63tq = stats.t.ppf(0.975, dfw)64lo, hi = eff - tq * se, eff + tq * se6566# (iii) OOS67wmeans = train.groupby("wd")["dy"].mean()68pred = test["wd"].map(wmeans)69rmse0 = float(np.sqrt(np.mean(test["dy"] ** 2)))70rmse1 = float(np.sqrt(np.mean((test["dy"] - pred) ** 2)))7172if p >= 0.05: concl = "no_detectable_effect"73elif lo > 0: concl = "positive_effect"74elif hi < 0: concl = "negative_effect"75else: concl = "no_detectable_effect"7677# figure78fig, ax = plt.subplots(figsize=(6, 4))79ax.errorbar(range(5), tab.mean_bp, yerr=1.96 * tab.se_bp, fmt="o", capsize=5)80ax.axhline(0, color="grey", lw=0.8)81ax.set_xticks(range(5)); ax.set_xticklabels([d[:3] for d in DAYS])82ax.set_ylabel("Mean daily change in DGS10 (bp)")83ax.set_title(f"10y Treasury yield, 2000-2019 (n={len(ch)}); 95% CI")84fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120); plt.close(fig)858687# ---- robustness / diagnostics (reported, do not change the primary conclusion) ----88import statsmodels.formula.api as smf89ch["mon"] = (ch.wd == "Monday").astype(int)90hc = smf.ols("dy ~ mon", ch).fit(cov_type="HC1")91pooled = smf.ols("dy ~ mon", ch).fit()92prev_date = df["date"].shift().iloc[1:]93ch["gap_days"] = (ch["date"] - prev_date).dt.days94noh = ch[ch.gap_days <= 3]95def anova_p(s): return float(stats.f_oneway(*[s.loc[s.wd == d, "dy"].values for d in DAYS]).pvalue)96imin = ch["dy"].idxmin(); imax = ch["dy"].idxmax()97robust = {98 "robust_mon_hc1_ci": [float(x) for x in hc.conf_int().loc["mon"]],99 "robust_mon_pooled_ci": [float(x) for x in pooled.conf_int().loc["mon"]],100 "robust_kruskal_p": float(stats.kruskal(*groups).pvalue),101 "diag_levene_p": float(stats.levene(*groups).pvalue),102 "diag_sd_by_weekday_bp": {d: float(g.std(ddof=1)) for d, g in zip(DAYS, groups)},103 "diag_lag1_autocorr": float(ch["dy"].autocorr()),104 "robust_anova_p_excl_post_holiday": anova_p(noh),105 "robust_n_excl_post_holiday": int(len(noh)),106 "robust_anova_p_train": anova_p(train), "robust_anova_p_test": anova_p(test),107 "diag_dy_sd_bp": float(ch["dy"].std(ddof=1)),108 "diag_dy_min_bp": [float(ch.loc[imin, "dy"]), str(ch.loc[imin, "date"].date())],109 "diag_dy_max_bp": [float(ch.loc[imax, "dy"]), str(ch.loc[imax, "date"].date())],110}111112res = {113 "task_id": "F6-conf-01",114 "n_missing_rows": n_missing, "n_train": int(len(train)), "n_test": int(len(test)),115 "anova_f_stat": float(F), "anova_p_value": float(p),116 "effect_estimate": float(eff), "effect_ci_low": float(lo), "effect_ci_high": float(hi),117 "effect_se": float(se), "effect_welch_df": float(dfw),118 "baseline_metric": rmse0, "model_metric": rmse1,119 "model_minus_baseline_rmse": rmse1 - rmse0,120 "train_weekday_means_bp": {d: float(wmeans[d]) for d in DAYS},121 "conclusion": concl, "feasibility": "complete", "feasibility_reason": None,122 "data_issues": data_issues, "n_input_rows": {"DGS10.csv": n_raw},123 "n_analysed": int(len(ch)), "random_seed": SEED,124 "input_sha256": {"DGS10.csv": sha},125 "first_change_date": str(ch.date.min().date()), "last_change_date": str(ch.date.max().date()),126 "last_train_date": str(train.date.max().date()), "first_test_date": str(test.date.min().date()),127}128res.update(robust)129with open(f"{OUT}/results.json", "w") as f: json.dump(res, f, indent=2)130print(json.dumps({k: res[k] for k in ["anova_f_stat", "anova_p_value", "effect_estimate", "effect_ci_low",131 "effect_ci_high", "baseline_metric", "model_metric", "conclusion"]}, indent=1))132