1"""F5-conf-02: Finite-sample bias of OLS AR(1) coefficient. Reads only data/design.csv."""2import json, hashlib, os3import numpy as np4import pandas as pd56ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))7DATA = os.path.join(ROOT, "data", "design.csv")8OUT = os.path.join(ROOT, "out")910raw = open(DATA, "rb").read()11sha = hashlib.sha256(raw).hexdigest()12d = pd.read_csv(DATA)13n_rows = int(len(d))14issues = []15if d["parameter"].duplicated().any():16 issues.append({"issue": "duplicated parameter names", "handling": "first occurrence used"})17p = dict(zip(d["parameter"].str.strip(), d["value"].astype(str).str.strip()))18required = ["n", "rho", "sigma", "intercept", "initial_condition", "n_replications", "seed"]19missing = [k for k in required if k not in p]20if missing:21 raise SystemExit(f"design missing parameters: {missing}")2223n = int(p["n"]); rho = float(p["rho"]); sigma = float(p["sigma"])24R = int(p["n_replications"]); seed = int(p["seed"])25intercept = p["intercept"].lower() == "yes"26assert intercept, "design requires intercept"27assert p["initial_condition"].lower() == "stationary"28assert abs(rho) < 1 and sigma > 0 and n >= 3 and R >= 22930# ---- analytical approximation and baseline ----31bias_analytical = -(1 + 3 * rho) / n32bias_analytical_n_minus_1 = -(1 + 3 * rho) / (n - 1) # sensitivity: effective sample size33se_asymptotic = np.sqrt((1 - rho ** 2) / n)3435# ---- Monte Carlo (vectorised, c = 0) ----36rng = np.random.default_rng(seed)37y = np.empty((R, n))38y[:, 0] = rng.normal(0.0, sigma / np.sqrt(1 - rho ** 2), size=R)39eps = rng.normal(0.0, sigma, size=(R, n - 1))40for t in range(1, n):41 y[:, t] = rho * y[:, t - 1] + eps[:, t - 1]42x = y[:, :-1]; z = y[:, 1:] # regress y_t on (1, y_{t-1}), t=2..n43xc = x - x.mean(axis=1, keepdims=True)44zc = z - z.mean(axis=1, keepdims=True)45rho_hat = (xc * zc).sum(axis=1) / (xc ** 2).sum(axis=1) # OLS slope with intercept4647finite = np.isfinite(rho_hat)48n_excl = int((~finite).sum())49rho_hat = rho_hat[finite]50n_an = int(rho_hat.size)51mean_rh = float(rho_hat.mean())52sd_rh = float(rho_hat.std(ddof=1))53bias_mc = mean_rh - rho54bias_mc_se = float(sd_rh / np.sqrt(n_an))55disc = float((bias_mc - bias_analytical) / bias_mc_se)5657# extra diagnostics (not required keys; written to separate file)58diag = {59 "design_sha256": sha,60 "bias_analytical_n_minus_1": bias_analytical_n_minus_1,61 "discrepancy_se_units_n_minus_1": float((bias_mc - bias_analytical_n_minus_1) / bias_mc_se),62 "bias_mc_ci95": [bias_mc - 1.96 * bias_mc_se, bias_mc + 1.96 * bias_mc_se],63 "rho_hat_median": float(np.median(rho_hat)),64 "share_rho_hat_below_rho": float((rho_hat < rho).mean()),65 "ratio_sd_mc_to_se_asymptotic": float(sd_rh / se_asymptotic),66 "relative_error_analytical_vs_mc": float((bias_analytical - bias_mc) / bias_mc),67 "n_nonfinite_excluded": n_excl,68}6970results = {71 "task_id": "F5-conf-02",72 "bias_analytical": bias_analytical,73 "bias_mc": bias_mc,74 "bias_mc_se": bias_mc_se,75 "rho_hat_mean": mean_rh,76 "rho_hat_sd_mc": sd_rh,77 "discrepancy_se_units": disc,78 "se_asymptotic": float(se_asymptotic),79 "n_replications": R,80 "feasibility": "complete",81 "feasibility_reason": None,82 "data_issues": issues,83 "n_input_rows": {"design.csv": n_rows},84 "n_analysed": n_an,85 "random_seed": seed,86}87os.makedirs(OUT, exist_ok=True)88with open(os.path.join(OUT, "results.json"), "w") as f:89 json.dump(results, f, indent=2)90with open(os.path.join(OUT, "diagnostics.json"), "w") as f:91 json.dump(diag, f, indent=2)92pd.DataFrame([{"step": 1, "rule": "drop replications with non-finite rho_hat",93 "n_before": R, "n_after": n_an, "n_removed": n_excl,94 "rationale": "degenerate regressor (zero variance) would make OLS undefined"}]95 ).to_csv(os.path.join(OUT, "exclusions.csv"), index=False)96print(json.dumps(results, indent=2)); print(json.dumps(diag, indent=2))97