1"""F5-conf-02: Finite-sample bias of OLS AR(1) coefficient (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.4"""5import json, math, csv6from pathlib import Path7import numpy as np8import pandas as pd910ROOT = Path(__file__).resolve().parent.parent11d = pd.read_csv(ROOT / "data" / "design.csv")12p = dict(zip(d["parameter"], d["value"]))13n = int(p["n"]); rho = float(p["rho"]); sigma = float(p["sigma"])14R = int(p["n_replications"]); seed = int(p["seed"])15assert str(p["intercept"]).strip().lower() == "yes"16assert str(p["initial_condition"]).strip().lower() == "stationary"17assert abs(rho) < 11819# --- Baseline / analytical quantities ---20bias_analytical = -(1 + 3 * rho) / n21se_asymptotic = math.sqrt((1 - rho ** 2) / n)2223# --- Monte Carlo (vectorised over replications), c = 0 ---24rng = np.random.default_rng(seed)25# Draw order: first y_1 for all R replications, then the (R, n-1) innovation matrix.26y = np.empty((R, n))27y[:, 0] = rng.normal(0.0, sigma / math.sqrt(1 - rho ** 2), size=R)28eps = rng.normal(0.0, sigma, size=(R, n - 1))29for t in range(1, n):30 y[:, t] = rho * y[:, t - 1] + eps[:, t - 1]3132# OLS of y_t on (1, y_{t-1}), t = 2..n: slope = cov(x,z)/var(x) with demeaning33x = y[:, :-1]; z = y[:, 1:]34xc = x - x.mean(axis=1, keepdims=True)35zc = z - z.mean(axis=1, keepdims=True)36rho_hat = (xc * zc).sum(axis=1) / (xc ** 2).sum(axis=1)3738# Spot check against explicit least squares for the first replication39X0 = np.column_stack([np.ones(n - 1), x[0]])40b0 = np.linalg.lstsq(X0, z[0], rcond=None)[0]41assert abs(b0[1] - rho_hat[0]) < 1e-104243# Exclusion ledger: non-finite estimates (none expected)44finite = np.isfinite(rho_hat)45rh = rho_hat[finite]46with open(ROOT / "out" / "exclusions.csv", "w", newline="") as f:47 w = csv.writer(f)48 w.writerow(["step", "rule", "n_before", "n_after", "n_removed", "rationale"])49 w.writerow([1, "drop replications with non-finite rho_hat", R, int(finite.sum()),50 int(R - finite.sum()), "OLS slope undefined if regressor has zero variance (probability zero)"])5152rho_hat_mean = float(rh.mean())53rho_hat_sd = float(rh.std(ddof=1))54bias_mc = rho_hat_mean - rho55bias_mc_se = rho_hat_sd / math.sqrt(rh.size)56disc = (bias_mc - bias_analytical) / bias_mc_se5758res = {59 "task_id": "F5-conf-02",60 "bias_analytical": bias_analytical,61 "bias_mc": bias_mc,62 "bias_mc_se": bias_mc_se,63 "rho_hat_mean": rho_hat_mean,64 "rho_hat_sd_mc": rho_hat_sd,65 "discrepancy_se_units": disc,66 "se_asymptotic": se_asymptotic,67 "n_replications": R,68 "feasibility": "complete",69 "feasibility_reason": None,70 "data_issues": [],71 "n_input_rows": {"design.csv": int(len(d))},72 "n_analysed": int(rh.size),73 "random_seed": seed,74 # supplementary (not required) descriptive quantities75 "supp_bias_analytical_n_minus_1": -(1 + 3 * rho) / (n - 1),76 "supp_rho_hat_median": float(np.median(rh)),77 "supp_ratio_bias_mc_to_analytical": bias_mc / bias_analytical,78 "supp_ratio_sd_mc_to_se_asymptotic": rho_hat_sd / se_asymptotic,79}80with open(ROOT / "out" / "results.json", "w") as f:81 json.dump(res, f, indent=2)82for k, v in res.items():83 print(f"{k}: {v}")84