1"""F6-conf-01 replication entry point. Run from workspace root:2 ~/wp29/venv/bin/python out/analysis.py3Reads data/DGS10.csv only; writes out/results.json, out/table1.csv, out/figure1.png, out/exclusions.csv.4No randomness is used (random_seed = null)."""5import json, hashlib, os6import numpy as np, pandas as pd7from scipy import stats8import matplotlib9matplotlib.use("Agg")10import matplotlib.pyplot as plt1112DATA = "data/DGS10.csv"; OUT = "out"13os.makedirs(OUT, exist_ok=True)14DAYS = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]1516sha = hashlib.sha256(open(DATA, "rb").read()).hexdigest()17raw = pd.read_csv(DATA, dtype={"observation_date": str, "DGS10": 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"] = pd.to_numeric(d["DGS10"], errors="raise")24d = d.sort_values("date").reset_index(drop=True)25n_avail = len(d)26d["dy"] = 100.0 * d["y"].diff()27d["dy"] = d["dy"].round(6) # remove float noise from 2-decimal percent inputs28ch = d.iloc[1:].copy() # first available observation has no change29ch["wd"] = ch["date"].dt.weekday30assert ch["wd"].between(0, 4).all()31ch["weekday"] = ch["wd"].map(dict(enumerate(DAYS)))3233pd.DataFrame([34 ["1", "drop rows with empty DGS10 (US market holidays)", n_raw, n_avail, n_raw - n_avail,35 "no yield published; required by TASK information rules"],36 ["2", "first available observation has no change", n_avail, len(ch), n_avail - len(ch),37 "change undefined for first observation"],38], columns=["step", "rule", "n_before", "n_after", "n_removed", "rationale"]).to_csv(39 f"{OUT}/exclusions.csv", index=False)4041train = ch[ch["date"] <= "2009-12-31"]42test = ch[ch["date"] >= "2010-01-01"]43assert len(train) + len(test) == len(ch)4445# (i) classical one-way ANOVA, all changes46groups = [ch.loc[ch["wd"] == k, "dy"].values for k in range(5)]47F, p = stats.f_oneway(*groups)4849# (ii) Monday contrast, Welch two-sample t interval50mon = groups[0]; oth = ch.loc[ch["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)54df_w = (v1 + v2) ** 2 / (v1 ** 2 / (len(mon) - 1) + v2 ** 2 / (len(oth) - 1))55tq = stats.t.ppf(0.975, df_w)56lo, hi = est - tq * se, est + tq * se57# robustness (reported in stdout only, not a required key): pooled t and HC1 regression58import statsmodels.api as sm59X = sm.add_constant((ch["wd"] == 0).astype(float).values)60hc = sm.OLS(ch["dy"].values, X).fit(cov_type="HC1")61print(f"[robustness] HC1 regression Monday coef {hc.params[1]:.4f}, CI {hc.conf_int()[1]}")6263# (iii) out-of-sample64wd_means_train = train.groupby("wd")["dy"].mean()65fc = test["wd"].map(wd_means_train)66assert fc.notna().all()67rmse_base = float(np.sqrt(np.mean(test["dy"].values ** 2)))68rmse_model = float(np.sqrt(np.mean((test["dy"].values - fc.values) ** 2)))69rmse_trainmean = float(np.sqrt(np.mean((test["dy"].values - train["dy"].mean()) ** 2)))70print(f"[context] RMSE of constant estimation-period mean forecast: {rmse_trainmean:.4f}")71print(f"[context] estimation-period weekday means (bp): {wd_means_train.round(4).to_dict()}")7273# table + figure74tab = pd.DataFrame({75 "weekday": DAYS,76 "n": [len(g) for g in groups],77 "mean_bp": [g.mean() for g in groups],78 "se_bp": [g.std(ddof=1) / np.sqrt(len(g)) for g in groups],79})80tab.to_csv(f"{OUT}/table1.csv", index=False)81fig, ax = plt.subplots(figsize=(6, 4))82ci = [stats.t.ppf(0.975, n - 1) * s for n, s in zip(tab.n, tab.se_bp)]83ax.errorbar(range(5), tab.mean_bp, yerr=ci, fmt="o", capsize=5)84ax.axhline(0, color="grey", lw=0.8)85ax.set_xticks(range(5)); ax.set_xticklabels(DAYS)86ax.set_ylabel("Mean one-day change in DGS10 (bp)")87ax.set_title(f"10y UST daily change by weekday, 2000-2019 (n={len(ch)})\n95% CIs; ANOVA F={F:.2f}, p={p:.3f}")88fig.tight_layout(); fig.savefig(f"{OUT}/figure1.png", dpi=120); plt.close(fig)8990if p < 0.05:91 conclusion = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")92else:93 conclusion = "no_detectable_effect"9495gaps = d["date"].diff().dt.days.iloc[1:].values96normal = np.where(ch["wd"].values == 0, 3, 1)97n_span = int((gaps > normal).sum())98res = {99 "task_id": "F6-conf-01",100 "n_missing_rows": n_missing,101 "n_train": int(len(train)), "n_test": int(len(test)),102 "anova_f_stat": float(F), "anova_p_value": float(p),103 "effect_estimate": float(est), "effect_ci_low": float(lo), "effect_ci_high": float(hi),104 "baseline_metric": rmse_base, "model_metric": rmse_model,105 "conclusion": conclusion,106 "feasibility": "complete", "feasibility_reason": None,107 "data_issues": [108 {"issue": f"{n_missing} rows with empty DGS10 field (US market holidays, incl. 2001-09-11/12)",109 "handling": "dropped before differencing, per TASK rules"},110 {"issue": f"{n_span} changes span a holiday (gap longer than the normal 1 day, or 3 days for Mondays); Monday changes always span a weekend",111 "handling": "kept and attributed to the weekday of t, as specified; no rescaling by elapsed days"},112 {"issue": "large outliers, e.g. -51 bp on 2009-03-18 (FOMC QE1 announcement)",113 "handling": "kept; verified as genuine market moves, not data errors"},114 ],115 "n_input_rows": {"DGS10.csv": int(n_raw)},116 "n_analysed": int(len(ch)),117 "random_seed": None,118 "monday_contrast_method": "Welch two-sample t interval (unequal variances), Monday vs Tue-Fri, all changes",119 "monday_contrast_se": float(se), "monday_contrast_welch_df": float(df_w),120 "robust_hc1_ci_low": float(hc.conf_int()[1][0]), "robust_hc1_ci_high": float(hc.conf_int()[1][1]),121 "rmse_ratio_model_over_baseline": rmse_model / rmse_base,122 "rmse_constant_train_mean": rmse_trainmean,123 "data_sha256": sha,124}125with open(f"{OUT}/results.json", "w") as f:126 json.dump(res, f, indent=2)127print(json.dumps({k: res[k] for k in ["n_train", "n_test", "anova_f_stat", "anova_p_value", "effect_estimate",128 "effect_ci_low", "effect_ci_high", "baseline_metric", "model_metric", "conclusion"]}, indent=1))129