"""F6-conf-01: weekday effects in daily changes of the 10y Treasury yield (DGS10). Run from workspace root: ~/wp29/venv/bin/python out/analysis.py""" import json, hashlib, os import numpy as np, pandas as pd from scipy import stats import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt SEED = None # no randomness used DATA = "data/DGS10.csv"; OUT = "out" os.makedirs(OUT, exist_ok=True) DAYS = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"] raw = pd.read_csv(DATA, dtype=str, keep_default_na=False) n_raw = len(raw) sha = hashlib.sha256(open(DATA, "rb").read()).hexdigest() raw["date"] = pd.to_datetime(raw["observation_date"], format="%Y-%m-%d") empty = raw["DGS10"].str.strip() == "" n_missing = int(empty.sum()) df = raw.loc[~empty].copy() df["y"] = pd.to_numeric(df["DGS10"], errors="raise") df = df.sort_values("date").reset_index(drop=True) # audit facts data_issues = [] n_dup = int(raw["observation_date"].duplicated().sum()) n_weekend = int((raw["date"].dt.dayofweek >= 5).sum()) data_issues.append({"issue": f"{n_missing} rows with empty DGS10 field (US market holidays)", "handling": "dropped before differencing, per task rule; change after a holiday spans the gap"}) if n_dup: data_issues.append({"issue": f"{n_dup} duplicated dates", "handling": "reported only"}) if n_weekend: data_issues.append({"issue": f"{n_weekend} weekend dates", "handling": "reported only"}) df["dy"] = 100.0 * df["y"].diff() df["dy"] = df["dy"].round(10) # remove float noise from percent->bp conversion ch = df.iloc[1:].copy() # first available obs has no change ch["wd"] = ch["date"].dt.day_name() assert set(ch["wd"]) <= set(DAYS) train = ch[ch["date"] <= "2009-12-31"]; test = ch[ch["date"] >= "2010-01-01"] excl = pd.DataFrame([ ["1", "drop rows with empty DGS10", n_raw, n_raw - n_missing, n_missing, "no published yield (holiday); task rule"], ["2", "first available observation has no change", n_raw - n_missing, len(ch), 1, "differencing"], ], columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]) excl.to_csv(f"{OUT}/exclusions.csv", index=False) # (i) ANOVA, all changes groups = [ch.loc[ch.wd == d, "dy"].values for d in DAYS] F, p = stats.f_oneway(*groups) # table tab = pd.DataFrame({"weekday": DAYS, "n": [len(g) for g in groups], "mean_bp": [g.mean() for g in groups], "se_bp": [g.std(ddof=1) / np.sqrt(len(g)) for g in groups]}) tab.to_csv(f"{OUT}/table1.csv", index=False) # (ii) Monday contrast, Welch two-sample mon = ch.loc[ch.wd == "Monday", "dy"].values; oth = ch.loc[ch.wd != "Monday", "dy"].values eff = mon.mean() - oth.mean() v1, v2 = mon.var(ddof=1) / len(mon), oth.var(ddof=1) / len(oth) se = np.sqrt(v1 + v2) dfw = (v1 + v2) ** 2 / (v1 ** 2 / (len(mon) - 1) + v2 ** 2 / (len(oth) - 1)) tq = stats.t.ppf(0.975, dfw) lo, hi = eff - tq * se, eff + tq * se # (iii) OOS wmeans = train.groupby("wd")["dy"].mean() pred = test["wd"].map(wmeans) rmse0 = float(np.sqrt(np.mean(test["dy"] ** 2))) rmse1 = float(np.sqrt(np.mean((test["dy"] - pred) ** 2))) if p >= 0.05: concl = "no_detectable_effect" elif lo > 0: concl = "positive_effect" elif hi < 0: concl = "negative_effect" else: concl = "no_detectable_effect" # figure fig, ax = plt.subplots(figsize=(6, 4)) ax.errorbar(range(5), tab.mean_bp, yerr=1.96 * tab.se_bp, fmt="o", capsize=5) ax.axhline(0, color="grey", lw=0.8) ax.set_xticks(range(5)); ax.set_xticklabels([d[:3] for d in DAYS]) ax.set_ylabel("Mean daily change in DGS10 (bp)") ax.set_title(f"10y Treasury yield, 2000-2019 (n={len(ch)}); 95% CI") fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120); plt.close(fig) # ---- robustness / diagnostics (reported, do not change the primary conclusion) ---- import statsmodels.formula.api as smf ch["mon"] = (ch.wd == "Monday").astype(int) hc = smf.ols("dy ~ mon", ch).fit(cov_type="HC1") pooled = smf.ols("dy ~ mon", ch).fit() prev_date = df["date"].shift().iloc[1:] ch["gap_days"] = (ch["date"] - prev_date).dt.days noh = ch[ch.gap_days <= 3] def anova_p(s): return float(stats.f_oneway(*[s.loc[s.wd == d, "dy"].values for d in DAYS]).pvalue) imin = ch["dy"].idxmin(); imax = ch["dy"].idxmax() robust = { "robust_mon_hc1_ci": [float(x) for x in hc.conf_int().loc["mon"]], "robust_mon_pooled_ci": [float(x) for x in pooled.conf_int().loc["mon"]], "robust_kruskal_p": float(stats.kruskal(*groups).pvalue), "diag_levene_p": float(stats.levene(*groups).pvalue), "diag_sd_by_weekday_bp": {d: float(g.std(ddof=1)) for d, g in zip(DAYS, groups)}, "diag_lag1_autocorr": float(ch["dy"].autocorr()), "robust_anova_p_excl_post_holiday": anova_p(noh), "robust_n_excl_post_holiday": int(len(noh)), "robust_anova_p_train": anova_p(train), "robust_anova_p_test": anova_p(test), "diag_dy_sd_bp": float(ch["dy"].std(ddof=1)), "diag_dy_min_bp": [float(ch.loc[imin, "dy"]), str(ch.loc[imin, "date"].date())], "diag_dy_max_bp": [float(ch.loc[imax, "dy"]), str(ch.loc[imax, "date"].date())], } res = { "task_id": "F6-conf-01", "n_missing_rows": n_missing, "n_train": int(len(train)), "n_test": int(len(test)), "anova_f_stat": float(F), "anova_p_value": float(p), "effect_estimate": float(eff), "effect_ci_low": float(lo), "effect_ci_high": float(hi), "effect_se": float(se), "effect_welch_df": float(dfw), "baseline_metric": rmse0, "model_metric": rmse1, "model_minus_baseline_rmse": rmse1 - rmse0, "train_weekday_means_bp": {d: float(wmeans[d]) for d in DAYS}, "conclusion": concl, "feasibility": "complete", "feasibility_reason": None, "data_issues": data_issues, "n_input_rows": {"DGS10.csv": n_raw}, "n_analysed": int(len(ch)), "random_seed": SEED, "input_sha256": {"DGS10.csv": sha}, "first_change_date": str(ch.date.min().date()), "last_change_date": str(ch.date.max().date()), "last_train_date": str(train.date.max().date()), "first_test_date": str(test.date.min().date()), } res.update(robust) with open(f"{OUT}/results.json", "w") as f: json.dump(res, f, indent=2) print(json.dumps({k: res[k] for k in ["anova_f_stat", "anova_p_value", "effect_estimate", "effect_ci_low", "effect_ci_high", "baseline_metric", "model_metric", "conclusion"]}, indent=1))