1"""Replication entry point for F6-conf-01: weekday effects in daily DGS10 changes."""2import json, os3import numpy as np, pandas as pd4from scipy import stats5import statsmodels.api as sm6import matplotlib7matplotlib.use("Agg")8import matplotlib.pyplot as plt910ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))11OUT = os.path.join(ROOT, "out")12raw = pd.read_csv(os.path.join(ROOT, "data", "DGS10.csv"), dtype=str, keep_default_na=False)13n_input = len(raw)14issues = []15raw["date"] = pd.to_datetime(raw["observation_date"])16miss = raw["DGS10"].str.strip() == ""17n_missing = int(miss.sum())18issues.append({"issue": f"{n_missing} rows with empty DGS10 field (US market holidays)",19 "handling": "dropped before differencing; change computed between consecutive available observations (so post-holiday changes span >1 business day)"})20df = raw.loc[~miss].copy()21df["y"] = pd.to_numeric(df["DGS10"], errors="coerce")22bad = int(df["y"].isna().sum())23if bad:24 issues.append({"issue": f"{bad} non-numeric DGS10 values", "handling": "dropped"})25 df = df.dropna(subset=["y"])26dup = int(df["date"].duplicated().sum())27if dup:28 issues.append({"issue": f"{dup} duplicated dates", "handling": "kept first"})29 df = df.drop_duplicates("date")30if not df["date"].is_monotonic_increasing:31 issues.append({"issue": "dates not sorted", "handling": "sorted"})32df = df.sort_values("date").reset_index(drop=True)33wkend = int((df["date"].dt.dayofweek >= 5).sum())34assert wkend == 0, "weekend observations present"3536df["dy"] = 100 * df["y"].diff()37df["dy"] = df["dy"].round(6) # remove float noise from decimal differencing38ch = df.dropna(subset=["dy"]).copy()39ch["wd"] = ch["date"].dt.dayofweek40names = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday"]41gap = (ch["date"] - df["date"].shift(1).loc[ch.index]).dt.days42n_long = int(((gap > 1) & ~((ch["wd"] == 0) & (gap == 3))).sum())43issues.append({"issue": f"{n_long} changes span a holiday gap (longer than the normal 1-day or Monday 3-day gap)",44 "handling": "kept, per task definition (attributed to weekday of t); noted as limitation"})4546train = ch[ch["date"] <= "2009-12-31"]; test = ch[ch["date"] >= "2010-01-01"]4748groups = [ch.loc[ch["wd"] == k, "dy"].values for k in range(5)]49F, p = stats.f_oneway(*groups)5051# Monday contrast: OLS on Monday indicator with HC3 robust SE52X = sm.add_constant((ch["wd"] == 0).astype(float).values)53fit = sm.OLS(ch["dy"].values, X).fit(cov_type="HC3")54eff = float(fit.params[1]); lo, hi = map(float, fit.conf_int(0.05)[1])55# Welch cross-check56w = stats.ttest_ind(groups[0], np.concatenate(groups[1:]), equal_var=False)5758tm = train.groupby("wd")["dy"].mean()59pred = test["wd"].map(tm)60rmse_base = float(np.sqrt(np.mean(test["dy"] ** 2)))61rmse_model = float(np.sqrt(np.mean((test["dy"] - pred) ** 2)))6263if p < 0.05:64 concl = "positive_effect" if lo > 0 else ("negative_effect" if hi < 0 else "no_detectable_effect")65else:66 concl = "no_detectable_effect"6768tab = pd.DataFrame({"weekday": names,69 "n": [len(g) for g in groups],70 "mean_bp": [g.mean() for g in groups],71 "se_bp": [g.std(ddof=1) / np.sqrt(len(g)) for g in groups]})72tab.to_csv(os.path.join(OUT, "table1.csv"), index=False, float_format="%.4f")7374crit = stats.t.ppf(0.975, tab["n"] - 1)75fig, ax = plt.subplots(figsize=(6, 4))76ax.errorbar(range(5), tab["mean_bp"], yerr=crit * tab["se_bp"], fmt="o", capsize=5)77ax.axhline(0, color="grey", lw=0.8)78ax.set_xticks(range(5)); ax.set_xticklabels(names)79ax.set_ylabel("Mean daily change in 10y yield (bp)")80ax.set_title(f"DGS10 daily change by weekday, 2000-2019 (95% CI)\nANOVA F={F:.2f}, p={p:.3f}")81fig.tight_layout(); fig.savefig(os.path.join(OUT, "figure1.png"), dpi=150)8283res = {"task_id": "F6-conf-01", "n_missing_rows": n_missing,84 "n_train": int(len(train)), "n_test": int(len(test)),85 "anova_f_stat": float(F), "anova_p_value": float(p),86 "effect_estimate": eff, "effect_ci_low": lo, "effect_ci_high": hi,87 "baseline_metric": rmse_base, "model_metric": rmse_model,88 "conclusion": concl, "feasibility": "complete", "feasibility_reason": None,89 "data_issues": issues, "n_input_rows": {"DGS10.csv": n_input},90 "n_analysed": int(len(ch)), "random_seed": None,91 "supplementary": {"welch_t": float(w.statistic), "welch_p": float(w.pvalue),92 "monday_contrast_hc3_p": float(fit.pvalues[1]),93 "train_weekday_means_bp": {names[k]: float(v) for k, v in tm.items()},94 "rmse_ratio_model_over_baseline": rmse_model / rmse_base,95 "test_mean_change_bp": float(test["dy"].mean())}}96with open(os.path.join(OUT, "results.json"), "w") as f:97 json.dump(res, f, indent=2)98print(json.dumps(res, indent=2)); print(tab)99