1"""F5-conf-02: finite-sample bias of OLS rho_hat in a Gaussian AR(1) with intercept.2Run from workspace root: ~/wp29/venv/bin/python out/analysis.py3Reads only data/design.csv; writes out/results.json and out/exclusions.csv deterministically."""4import json, hashlib, os5import numpy as np6import pandas as pd78ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))9DESIGN = os.path.join(ROOT, "data", "design.csv")10OUT = os.path.join(ROOT, "out")1112d = pd.read_csv(DESIGN)13p = dict(zip(d["parameter"].astype(str).str.strip(), d["value"].astype(str).str.strip()))14issues = []15n = int(p["n"]); rho = float(p["rho"]); sigma = float(p["sigma"])16R = int(p["n_replications"]); seed = int(p["seed"])17intercept = p["intercept"].lower() == "yes"18init = p["initial_condition"].lower()19assert intercept, "design requires intercept"20assert init == "stationary" and abs(rho) < 1, "stationary init requires |rho|<1"21assert n >= 3 and R >= 2 and sigma > 02223# ---- P2 baseline quantities (closed form) ----24bias_analytical = -(1 + 3 * rho) / n25se_asymptotic = float(np.sqrt((1 - rho ** 2) / n))26# sensitivity: same formula with effective number of regression observations n-127bias_analytical_neff = -(1 + 3 * rho) / (n - 1)2829# ---- P3 Monte Carlo ----30# Draw order: one (R, n) standard-normal matrix; column 0 scaled to the stationary31# sd for y_1, columns 1..n-1 scaled by sigma as innovations e_2..e_n. c = 0.32rng = np.random.default_rng(seed)33z = rng.standard_normal((R, n))34y = np.empty((R, n))35y[:, 0] = z[:, 0] * sigma / np.sqrt(1 - rho ** 2)36eps = sigma * z[:, 1:]37for t in range(1, n):38 y[:, t] = rho * y[:, t - 1] + eps[:, t - 1]39x = y[:, :-1]; yy = y[:, 1:] # t = 2..n (n-1 pairs)40xc = x - x.mean(axis=1, keepdims=True)41yc = yy - yy.mean(axis=1, keepdims=True)42rho_hat = (xc * yc).sum(axis=1) / (xc ** 2).sum(axis=1)4344finite = np.isfinite(rho_hat)45n_bad = int((~finite).sum())46if n_bad:47 issues.append({"issue": f"{n_bad} non-finite rho_hat", "handling": "excluded; see exclusions.csv"})48rh = rho_hat[finite]49Rused = int(rh.size)5051rho_hat_mean = float(rh.mean())52rho_hat_sd = float(rh.std(ddof=1))53bias_mc = rho_hat_mean - rho54bias_mc_se = rho_hat_sd / np.sqrt(Rused)55disc = (bias_mc - bias_analytical) / bias_mc_se56disc_neff = (bias_mc - bias_analytical_neff) / bias_mc_se57q = np.quantile(rh, [0.025, 0.5, 0.975])5859pd.DataFrame([{"step": 1, "rule": "drop replications with non-finite rho_hat",60 "n_before": R, "n_after": Rused, "n_removed": R - Rused,61 "rationale": "OLS undefined if regressor has zero variance"}]62 ).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)6364res = {65 "task_id": "F5-conf-02",66 "bias_analytical": float(bias_analytical),67 "bias_mc": float(bias_mc),68 "bias_mc_se": float(bias_mc_se),69 "rho_hat_mean": rho_hat_mean,70 "rho_hat_sd_mc": rho_hat_sd,71 "discrepancy_se_units": float(disc),72 "se_asymptotic": se_asymptotic,73 "n_replications": R,74 "feasibility": "complete",75 "feasibility_reason": None,76 "data_issues": issues,77 "n_input_rows": {"design.csv": int(len(d))},78 "n_analysed": Rused,79 "random_seed": seed,80 # supplementary keys81 "design_n": n, "design_rho": rho, "design_sigma": sigma,82 "bias_mc_ci95": [float(bias_mc - 1.96 * bias_mc_se), float(bias_mc + 1.96 * bias_mc_se)],83 "bias_analytical_neff": float(bias_analytical_neff),84 "discrepancy_se_units_neff": float(disc_neff),85 "bias_ratio_mc_to_analytical": float(bias_mc / bias_analytical),86 "rho_hat_median": float(q[1]),87 "rho_hat_q025": float(q[0]), "rho_hat_q975": float(q[2]),88 "sd_ratio_mc_to_asymptotic": float(rho_hat_sd / se_asymptotic),89 "share_rho_hat_below_rho": float((rh < rho).mean()),90 "design_sha256": hashlib.sha256(open(DESIGN, "rb").read()).hexdigest(),91}92with open(os.path.join(OUT, "results.json"), "w") as f:93 json.dump(res, f, indent=2)94print(json.dumps(res, indent=2))95