1"""Replication script for F6-conf-01: day-of-week effects in daily 10y Treasury yield changes."""2import json, os3import numpy as np, pandas as pd4from scipy import stats5import matplotlib6matplotlib.use("Agg")7import matplotlib.pyplot as plt89OUT = "out"; os.makedirs(OUT, exist_ok=True)10raw = pd.read_csv("data/DGS10.csv", dtype=str, keep_default_na=False)11n_input = len(raw)12issues = []13raw["date"] = pd.to_datetime(raw["observation_date"], format="%Y-%m-%d")14empty = raw["DGS10"].str.strip() == ""15n_missing = int(empty.sum())16issues.append({"issue": f"{n_missing} rows with empty DGS10 field (US market holidays / no publication)",17 "handling": "rows dropped before computing changes, as required; the change after a gap spans >1 business day and is attributed to the weekday of its end date"})18vals = pd.to_numeric(raw.loc[~empty, "DGS10"], errors="coerce")19n_nonnum = int(vals.isna().sum())20if n_nonnum:21 issues.append({"issue": f"{n_nonnum} non-numeric DGS10 values", "handling": "treated as missing and dropped"})22df = raw.loc[~empty].assign(y=vals).dropna(subset=["y"]).sort_values("date")23n_dup = int(df["date"].duplicated().sum())24if n_dup:25 issues.append({"issue": f"{n_dup} duplicated dates", "handling": "kept first occurrence"})26 df = df.drop_duplicates("date")27n_weekend = int((df["date"].dt.dayofweek >= 5).sum())28if n_weekend:29 issues.append({"issue": f"{n_weekend} weekend observations", "handling": "dropped"})30 df = df[df["date"].dt.dayofweek < 5]31df["dy"] = 100 * df["y"].diff()32df = df.dropna(subset=["dy"]).copy()33gap = df["date"].diff().dt.days34df["wd"] = df["date"].dt.dayofweek35# count changes spanning a holiday (calendar gap larger than normal)36normal = np.where(df["wd"] == 0, 3, 1)37prev_dates = df["date"].shift(1)38span = (df["date"] - prev_dates).dt.days.values39long_span = int(np.nansum(span[1:] > normal[1:]))40issues.append({"issue": f"{long_span} changes span a holiday gap (longer than the usual 1 or 3 calendar days); this mixes horizons across weekdays (e.g. Tuesdays after Monday holidays)",41 "handling": "kept as specified by TASK.md; noted as a limitation"})42names = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]43df["weekday"] = df["wd"].map(dict(enumerate(names)))4445train = df[df["date"] <= "2009-12-31"]; test = df[df["date"] >= "2010-01-01"]4647groups = [df.loc[df["wd"] == k, "dy"].values for k in range(5)]48F, p = stats.f_oneway(*groups)4950mon = df.loc[df["wd"] == 0, "dy"].values; oth = df.loc[df["wd"] != 0, "dy"].values51est = mon.mean() - oth.mean()52v1, v2 = mon.var(ddof=1) / len(mon), oth.var(ddof=1) / len(oth)53se = np.sqrt(v1 + v2)54dof = (v1 + v2) ** 2 / (v1 ** 2 / (len(mon) - 1) + v2 ** 2 / (len(oth) - 1))55tq = stats.t.ppf(0.975, dof)56lo, hi = est - tq * se, est + tq * se57welch_p = float(2 * stats.t.sf(abs(est / se), dof))5859wmeans = train.groupby("wd")["dy"].mean()60pred = test["wd"].map(wmeans).values61rmse0 = float(np.sqrt(np.mean(test["dy"].values ** 2)))62rmse1 = float(np.sqrt(np.mean((test["dy"].values - pred) ** 2)))63# supplementary: training-mean (no weekday) forecast64rmse_mean = float(np.sqrt(np.mean((test["dy"].values - train["dy"].mean()) ** 2)))6566if p < 0.05:67 concl = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")68else:69 concl = "no_detectable_effect"7071tab = pd.DataFrame({"weekday": names, "n": [len(g) for g in groups],72 "mean_bp": [g.mean() for g in groups],73 "se_bp": [g.std(ddof=1) / np.sqrt(len(g)) for g in groups]})74tab.to_csv(f"{OUT}/table1.csv", index=False, float_format="%.6f")7576fig, ax = plt.subplots(figsize=(6, 4))77ci = [stats.t.ppf(0.975, n - 1) * s for n, s in zip(tab["n"], tab["se_bp"])]78ax.errorbar(range(5), tab["mean_bp"], yerr=ci, fmt="o", capsize=5)79ax.axhline(0, color="grey", lw=0.8, ls="--")80ax.set_xticks(range(5)); ax.set_xticklabels(names)81ax.set_ylabel("Mean one-day change in DGS10 (bp)")82ax.set_title(f"10y Treasury yield change by weekday, 2000-2019\n(95% CI; ANOVA F={F:.2f}, p={p:.3f})")83fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=150)8485res = {"task_id": "F6-conf-01", "n_missing_rows": n_missing,86 "n_train": int(len(train)), "n_test": int(len(test)),87 "anova_f_stat": float(F), "anova_p_value": float(p),88 "effect_estimate": float(est), "effect_ci_low": float(lo), "effect_ci_high": float(hi),89 "baseline_metric": rmse0, "model_metric": rmse1, "conclusion": concl,90 "feasibility": "complete", "feasibility_reason": None, "data_issues": issues,91 "n_input_rows": {"DGS10.csv": int(n_input)}, "n_analysed": int(len(df)), "random_seed": None,92 "supplementary": {"monday_contrast_method": "Welch two-sample t (Monday vs other weekdays)",93 "welch_se_bp": float(se), "welch_df": float(dof), "welch_p_value": welch_p,94 "train_weekday_means_bp": {names[k]: float(v) for k, v in wmeans.items()},95 "rmse_train_overall_mean_bp": rmse_mean,96 "rmse_ratio_model_to_baseline": rmse1 / rmse0,97 "n_valid_obs": int(len(df) + 1)}}98with open(f"{OUT}/results.json", "w") as f:99 json.dump(res, f, indent=2)100print(json.dumps(res, indent=2)); print(tab)101