SPB Git

spb/wp29_uqo

Public
6.0 KB · 125 lines python
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 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={"observation_date": str})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()19y = pd.to_numeric(raw["DGS10"], errors="coerce")20n_missing = int(y.isna().sum())21d = raw.loc[y.notna(), ["date"]].copy(); d["y"] = y[y.notna()].values22n_avail = len(d)23d["dy"] = 100.0 * d["y"].diff()            # bp, between consecutive available obs24d["gap_days"] = d["date"].diff().dt.days25ch = d.dropna(subset=["dy"]).copy()26ch["dy"] = ch["dy"].round(10)              # remove float noise from 2-decimal differences27ch["wd"] = ch["date"].dt.dayofweek28assert ch["wd"].between(0, 4).all()29n_changes = len(ch)3031pd.DataFrame([32    dict(step=1, rule="drop rows with empty DGS10 (holidays)", n_before=n_raw, n_after=n_avail,33         n_removed=n_raw - n_avail, rationale="no published yield; required by TASK"),34    dict(step=2, rule="first available observation has no change", n_before=n_avail, n_after=n_changes,35         n_removed=1, rationale="differencing"),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) == n_changes4041# (i) one-way ANOVA, all changes42groups = [ch.loc[ch.wd == k, "dy"].values for k in range(5)]43F, p = stats.f_oneway(*groups)4445# (ii) Monday contrast, Welch two-sample t interval46mon = ch.loc[ch.wd == 0, "dy"].values; oth = ch.loc[ch.wd != 0, "dy"].values47eff = mon.mean() - oth.mean()48v1, v2 = mon.var(ddof=1) / len(mon), oth.var(ddof=1) / len(oth)49se = np.sqrt(v1 + v2)50dfw = (v1 + v2) ** 2 / (v1 ** 2 / (len(mon) - 1) + v2 ** 2 / (len(oth) - 1))51tq = stats.t.ppf(0.975, dfw)52lo, hi = eff - tq * se, eff + tq * se53welch_p = 2 * stats.t.sf(abs(eff / se), dfw)5455# (iii) out-of-sample: weekday means from estimation period only56wd_means_train = train.groupby("wd")["dy"].mean()57pred = test["wd"].map(wd_means_train).values58rmse_rw = float(np.sqrt(np.mean(test["dy"].values ** 2)))59rmse_wd = float(np.sqrt(np.mean((test["dy"].values - pred) ** 2)))60# supplementary: Diebold-Mariano (squared loss, h=1, plain variance since h=1)61dl = test["dy"].values ** 2 - (test["dy"].values - pred) ** 2   # >0 means weekday model better62dm = dl.mean() / (dl.std(ddof=1) / np.sqrt(len(dl)))63dm_p = 2 * stats.norm.sf(abs(dm))64rmse_train_mean = float(np.sqrt(np.mean((test["dy"].values - train["dy"].mean()) ** 2)))6566# table + figure67tab = pd.DataFrame({"weekday": DAYS,68                    "n": [len(g) for g in groups],69                    "mean_bp": [g.mean() for g in groups],70                    "se_bp": [g.std(ddof=1) / np.sqrt(len(g)) for g in groups]})71tab.to_csv(f"{OUT}/table1.csv", index=False, float_format="%.6f")72crit = [stats.t.ppf(0.975, n - 1) for n in tab.n]73fig, ax = plt.subplots(figsize=(6, 4))74ax.errorbar(range(5), tab.mean_bp, yerr=np.array(crit) * tab.se_bp, fmt="o", capsize=5)75ax.axhline(0, color="grey", lw=0.8, ls="--")76ax.set_xticks(range(5)); ax.set_xticklabels(DAYS)77ax.set_ylabel("Mean daily change in DGS10 (bp)")78ax.set_title("10y Treasury yield: mean 1-day change by weekday\n2000-2019, 95% CI (t, per-weekday SE)")79fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=150); plt.close(fig)8081if p >= 0.05:82    concl = "no_detectable_effect"83elif lo > 0:84    concl = "positive_effect"85elif hi < 0:86    concl = "negative_effect"87else:88    concl = "no_detectable_effect"8990n_long_gap = int((ch["gap_days"] > np.where(ch["wd"] == 0, 3, 1)).sum())  # longer than a normal weekday/weekend gap91res = {92    "task_id": "F6-conf-01",93    "n_missing_rows": n_missing,94    "n_train": int(len(train)), "n_test": int(len(test)),95    "anova_f_stat": float(F), "anova_p_value": float(p),96    "effect_estimate": float(eff), "effect_ci_low": float(lo), "effect_ci_high": float(hi),97    "baseline_metric": rmse_rw, "model_metric": rmse_wd,98    "conclusion": concl,99    "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, as required by TASK"},103        {"issue": f"{n_long_gap} changes span a holiday (calendar gap longer than normal: > 3 days for Monday, > 1 day otherwise), e.g. Tuesday after a Monday holiday covers two trading sessions' worth of calendar time",104         "handling": "kept and attributed to the weekday of t per TASK rule; noted as a limitation"},105    ],106    "n_input_rows": {"DGS10.csv": int(n_raw)},107    "n_analysed": int(n_changes),108    "random_seed": None,109    "supplementary": {110        "effect_method": "Welch two-sample t interval (Monday vs other weekdays, all changes)",111        "effect_se": float(se), "effect_welch_df": float(dfw), "effect_welch_p_value": float(welch_p),112        "n_monday": int(len(mon)), "n_other": int(len(oth)),113        "train_weekday_means_bp": {DAYS[k]: float(v) for k, v in wd_means_train.items()},114        "rmse_ratio_model_over_baseline": rmse_wd / rmse_rw,115        "dm_stat_sqloss_positive_favours_weekday": float(dm), "dm_p_value": float(dm_p),116        "rmse_test_train_unconditional_mean": rmse_train_mean,117        "n_changes_spanning_holiday_gap": n_long_gap,118        "n_available_obs": int(n_avail),119    },120}121with open(f"{OUT}/results.json", "w") as f:122    json.dump(res, f, indent=2)123print(json.dumps({k: res[k] for k in ["n_train", "n_test", "anova_f_stat", "anova_p_value",124      "effect_estimate", "effect_ci_low", "effect_ci_high", "baseline_metric", "model_metric", "conclusion"]}, indent=1))125
Files · main
Menu

Theme