1"""F5-conf-02: finite-sample bias of OLS AR(1) coefficient (with intercept).2Reads only data/design.csv; writes out/results.json (and out/exclusions.csv). Deterministic."""3import json, os, hashlib4import numpy as np5import pandas as pd67ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))8DATA = os.path.join(ROOT, "data", "design.csv")9OUT = os.path.join(ROOT, "out")1011raw = pd.read_csv(DATA)12n_rows = int(len(raw))13d = dict(zip(raw["parameter"].astype(str).str.strip(), raw["value"].astype(str).str.strip()))14data_issues = []15if raw["parameter"].duplicated().any():16 data_issues.append({"issue": "duplicated parameter rows", "handling": "first occurrence used"})1718n = int(d["n"]); rho = float(d["rho"]); sigma = float(d["sigma"])19intercept = d["intercept"].lower() == "yes"20init = d["initial_condition"].lower()21R = int(d["n_replications"]); seed = int(d["seed"])22assert intercept and init == "stationary" and abs(rho) < 1 and n > 3 and sigma > 02324# --- baseline / analytical quantities ---25bias_analytical = -(1.0 + 3.0 * rho) / n26se_asymptotic = float(np.sqrt((1.0 - rho ** 2) / n))2728# --- Monte Carlo ---29# One (R, n) standard-normal draw: column 0 -> stationary y_1, columns 1..n-1 -> innovations.30rng = np.random.default_rng(seed)31z = rng.standard_normal((R, n))32y = np.empty((R, n))33y[:, 0] = sigma / np.sqrt(1.0 - rho ** 2) * z[:, 0]34eps = sigma * z[:, 1:]35for t in range(1, n):36 y[:, t] = rho * y[:, t - 1] + eps[:, t - 1] # c = 037x = y[:, :-1]; yy = y[:, 1:] # t = 2..n (n-1 obs)38xc = x - x.mean(axis=1, keepdims=True)39yc = yy - yy.mean(axis=1, keepdims=True)40sxx = (xc * xc).sum(axis=1)41rho_hat = (xc * yc).sum(axis=1) / sxx # OLS slope with intercept (FWL)4243ok = np.isfinite(rho_hat) & (sxx > 0)44n_bad = int((~ok).sum())45rows = [{"step": 1, "rule": "drop replications with non-finite rho_hat or zero regressor variance",46 "n_before": R, "n_after": R - n_bad, "n_removed": n_bad,47 "rationale": "OLS undefined; none expected under continuous Gaussian DGP"}]48pd.DataFrame(rows).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)49if n_bad:50 data_issues.append({"issue": f"{n_bad} degenerate replications", "handling": "excluded"})51rh = rho_hat[ok]52R_used = int(rh.size)5354rho_hat_mean = float(rh.mean())55rho_hat_sd = float(rh.std(ddof=1))56bias_mc = rho_hat_mean - rho57bias_mc_se = rho_hat_sd / np.sqrt(R_used)58disc = (bias_mc - bias_analytical) / bias_mc_se5960res = {61 "task_id": "F5-conf-02",62 "bias_analytical": float(bias_analytical),63 "bias_mc": float(bias_mc),64 "bias_mc_se": float(bias_mc_se),65 "rho_hat_mean": rho_hat_mean,66 "rho_hat_sd_mc": rho_hat_sd,67 "discrepancy_se_units": float(disc),68 "se_asymptotic": se_asymptotic,69 "n_replications": int(R),70 "feasibility": "complete",71 "feasibility_reason": None,72 "data_issues": data_issues,73 "n_input_rows": {"design.csv": n_rows},74 "n_analysed": R_used,75 "random_seed": seed,76 # supplementary diagnostics77 "supp_bias_analytical_neff": float(-(1.0 + 3.0 * rho) / (n - 1)),78 "supp_rho_hat_median": float(np.median(rh)),79 "supp_bias_mc_ci95_low": float(bias_mc - 1.96 * bias_mc_se),80 "supp_bias_mc_ci95_high": float(bias_mc + 1.96 * bias_mc_se),81 "supp_ratio_bias_mc_to_analytical": float(bias_mc / bias_analytical),82 "supp_ratio_sd_mc_to_se_asymptotic": float(rho_hat_sd / se_asymptotic),83 "supp_bias_mc_over_sd": float(bias_mc / rho_hat_sd),84 "supp_rmse_mc": float(np.sqrt(np.mean((rh - rho) ** 2))),85 "supp_share_rho_hat_below_rho": float(np.mean(rh < rho)),86 "design_sha256": hashlib.sha256(open(DATA, "rb").read()).hexdigest(),87}88with open(os.path.join(OUT, "results.json"), "w") as f:89 json.dump(res, f, indent=2)90print(json.dumps({k: res[k] for k in ["bias_analytical", "bias_mc", "bias_mc_se", "discrepancy_se_units",91 "rho_hat_mean", "rho_hat_sd_mc", "se_asymptotic", "n_analysed"]}, indent=2))92