1"""Replication entry point for F6-conf-01. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3Reads only data/DGS10.csv; writes out/results.json, out/table1.csv, out/figure1.png, out/exclusions.csv."""4import json, os5import numpy as np6import pandas as pd7from scipy import stats8import matplotlib9matplotlib.use("Agg")10import matplotlib.pyplot as plt1112OUT = "out"13os.makedirs(OUT, exist_ok=True)14SEED = None # no randomness used15DAYS = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]1617raw = pd.read_csv("data/DGS10.csv", dtype=str, keep_default_na=False)18n_raw = len(raw)19empty = raw["DGS10"].str.strip() == ""20n_missing = int(empty.sum())21d = raw.loc[~empty].copy()22d["date"] = pd.to_datetime(d["observation_date"])23d["y"] = d["DGS10"].astype(float)24assert d["date"].is_monotonic_increasing and not d["date"].duplicated().any()25d = d.reset_index(drop=True)26d["dy"] = 100.0 * d["y"].diff()27d["gap_days"] = d["date"].diff().dt.days28ch = d.iloc[1:].copy() # first available obs has no change29ch["wd"] = ch["date"].dt.dayofweek30assert ch["wd"].between(0, 4).all()31ch["weekday"] = ch["wd"].map(dict(enumerate(DAYS)))32n_ch = len(ch)3334excl = pd.DataFrame([35 {"step": 1, "rule": "drop rows with empty DGS10 (market holidays)", "n_before": n_raw,36 "n_after": n_raw - n_missing, "n_removed": n_missing, "rationale": "no yield published; required by TASK info rules"},37 {"step": 2, "rule": "first available observation has no change", "n_before": n_raw - n_missing,38 "n_after": n_ch, "n_removed": n_raw - n_missing - n_ch, "rationale": "difference undefined"},39])40excl.to_csv(os.path.join(OUT, "exclusions.csv"), index=False)4142# (i) ANOVA, all changes43groups = [ch.loc[ch["wd"] == k, "dy"].values for k in range(5)]44F, p = stats.f_oneway(*groups)4546# table 147tab = pd.DataFrame({"weekday": DAYS,48 "n": [len(g) for g in groups],49 "mean_bp": [g.mean() for g in groups],50 "se_bp": [g.std(ddof=1) / np.sqrt(len(g)) for g in groups]})51tab.to_csv(os.path.join(OUT, "table1.csv"), index=False)5253# (ii) Monday contrast, Welch two-sample interval54mon, oth = groups[0], np.concatenate(groups[1:])55eff = mon.mean() - oth.mean()56v1, v2 = mon.var(ddof=1) / len(mon), oth.var(ddof=1) / len(oth)57se = np.sqrt(v1 + v2)58df_w = (v1 + v2) ** 2 / (v1 ** 2 / (len(mon) - 1) + v2 ** 2 / (len(oth) - 1))59tcrit = stats.t.ppf(0.975, df_w)60lo, hi = eff - tcrit * se, eff + tcrit * se61welch_p = 2 * stats.t.sf(abs(eff / se), df_w)6263# (iii) out-of-sample64cut = pd.Timestamp("2009-12-31")65tr, te = ch[ch["date"] <= cut], ch[ch["date"] > cut]66wd_means = tr.groupby("wd")["dy"].mean()67assert tr["date"].max() < te["date"].min()68pred = te["wd"].map(wd_means).values69rmse0 = float(np.sqrt(np.mean(te["dy"].values ** 2)))70rmse1 = float(np.sqrt(np.mean((te["dy"].values - pred) ** 2)))71# supplementary: unconditional train-mean forecast72rmse_mean = float(np.sqrt(np.mean((te["dy"].values - tr["dy"].mean()) ** 2)))7374if p >= 0.05:75 concl = "no_detectable_effect"76elif lo > 0:77 concl = "positive_effect"78elif hi < 0:79 concl = "negative_effect"80else:81 concl = "no_detectable_effect"8283n_long_gap = int((ch["gap_days"] > 4).sum())84n_holiday_spanning = int(((ch["gap_days"] > 1) & ~((ch["wd"] == 0) & (ch["gap_days"] == 3))).sum())85issues = [86 {"issue": f"{n_missing} rows with empty DGS10 (US market holidays)",87 "handling": "dropped before differencing, as required; changes computed between consecutive available observations"},88 {"issue": f"{n_holiday_spanning} changes span a holiday (gap longer than a normal weekday/weekend gap); {n_long_gap} span more than 4 calendar days",89 "handling": "kept and attributed to the weekday of t per the task rules; noted as limitation (multi-day changes have larger variance)"},90]9192res = {93 "task_id": "F6-conf-01",94 "feasibility": "complete",95 "feasibility_reason": None,96 "data_issues": issues,97 "n_input_rows": {"DGS10.csv": int(n_raw)},98 "n_analysed": int(n_ch),99 "random_seed": SEED,100 "n_missing_rows": n_missing,101 "n_train": int(len(tr)),102 "n_test": int(len(te)),103 "anova_f_stat": float(F),104 "anova_p_value": float(p),105 "effect_estimate": float(eff),106 "effect_ci_low": float(lo),107 "effect_ci_high": float(hi),108 "baseline_metric": rmse0,109 "model_metric": rmse1,110 "conclusion": concl,111 # supplementary diagnostics112 "supp_effect_se": float(se),113 "supp_effect_welch_df": float(df_w),114 "supp_effect_welch_p": float(welch_p),115 "supp_rmse_ratio_model_over_baseline": rmse1 / rmse0,116 "supp_rmse_train_mean_forecast": rmse_mean,117 "supp_levene_p": float(stats.levene(*groups, center="median").pvalue),118 "supp_kruskal_p": float(stats.kruskal(*groups).pvalue),119 "supp_train_weekday_means_bp": {DAYS[k]: float(v) for k, v in wd_means.items()},120 "supp_n_holiday_spanning_changes": n_holiday_spanning,121}122with open(os.path.join(OUT, "results.json"), "w") as f:123 json.dump(res, f, indent=2)124125# figure 1126ci = stats.t.ppf(0.975, tab["n"] - 1) * tab["se_bp"]127fig, ax = plt.subplots(figsize=(6, 4))128ax.errorbar(range(5), tab["mean_bp"], yerr=ci, fmt="o", capsize=5)129ax.axhline(0, color="grey", lw=0.8)130ax.set_xticks(range(5)); ax.set_xticklabels(DAYS)131ax.set_ylabel("Mean one-day change in DGS10 (bp)")132ax.set_title("10y Treasury yield: mean daily change by weekday, 2000-2019\n(95% CI)")133fig.tight_layout(); fig.savefig(os.path.join(OUT, "figure1.png"), dpi=120)134135print(json.dumps({k: res[k] for k in ["n_analysed", "n_train", "n_test", "anova_f_stat", "anova_p_value",136 "effect_estimate", "effect_ci_low", "effect_ci_high", "baseline_metric", "model_metric", "conclusion"]}, indent=1))137print(tab.to_string(index=False))138