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."""4import json, csv, os5import numpy as np6import pandas as pd78ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))9d = pd.read_csv(os.path.join(ROOT, "data", "design.csv"))10p = dict(zip(d["parameter"], d["value"]))11n = int(p["n"]); rho = float(p["rho"]); sigma = float(p["sigma"])12R = int(p["n_replications"]); seed = int(p["seed"])13assert str(p["intercept"]).strip().lower() == "yes"14assert str(p["initial_condition"]).strip().lower() == "stationary"15assert abs(rho) < 11617data_issues = []18if d["parameter"].duplicated().any():19 data_issues.append({"issue": "duplicated parameter names", "handling": "first kept"})2021# --- Simulation (vectorised). RNG draw order: first R stationary initial values,22# then an R x (n-1) matrix of innovations. c = 0.23rng = np.random.default_rng(seed)24y = np.empty((R, n))25y[:, 0] = rng.standard_normal(R) * sigma / np.sqrt(1 - rho**2)26eps = rng.standard_normal((R, n - 1)) * sigma27for t in range(1, n):28 y[:, t] = rho * y[:, t - 1] + eps[:, t - 1]2930# --- OLS of y_t on (1, y_{t-1}), t = 2..n: slope = cov / var (demeaned)31x = y[:, :-1]; z = y[:, 1:]32xd = x - x.mean(axis=1, keepdims=True); zd = z - z.mean(axis=1, keepdims=True)33sxx = (xd**2).sum(axis=1)34ok = np.isfinite(sxx) & (sxx > 0)35rho_hat = (xd * zd).sum(axis=1)[ok] / sxx[ok]3637# cross-check a few replications against numpy lstsq38for i in range(5):39 X = np.column_stack([np.ones(n - 1), x[i]])40 b = np.linalg.lstsq(X, z[i], rcond=None)[0]41 assert abs(b[1] - rho_hat[i]) < 1e-104243with open(os.path.join(ROOT, "out", "exclusions.csv"), "w", newline="") as f:44 w = csv.writer(f)45 w.writerow(["step", "rule", "n_before", "n_after", "n_removed", "rationale"])46 w.writerow([1, "drop replications with zero/non-finite variance of y_{t-1}",47 R, int(ok.sum()), int(R - ok.sum()), "OLS slope undefined"])4849n_an = int(rho_hat.size)50mean = float(rho_hat.mean()); sd = float(rho_hat.std(ddof=1))51bias_mc = mean - rho; se = sd / np.sqrt(n_an)52bias_an = -(1 + 3 * rho) / n53res = {54 "task_id": "F5-conf-02",55 "bias_analytical": bias_an,56 "bias_mc": bias_mc,57 "bias_mc_se": se,58 "rho_hat_mean": mean,59 "rho_hat_sd_mc": sd,60 "discrepancy_se_units": (bias_mc - bias_an) / se,61 "se_asymptotic": float(np.sqrt((1 - rho**2) / n)),62 "n_replications": R,63 "feasibility": "complete",64 "feasibility_reason": None,65 "data_issues": data_issues,66 "n_input_rows": {"design.csv": int(len(d))},67 "n_analysed": n_an,68 "random_seed": seed,69 # supplementary (not required)70 "supp_bias_analytical_n_minus_1": -(1 + 3 * rho) / (n - 1),71 "supp_bias_mc_ci95": [bias_mc - 1.96 * se, bias_mc + 1.96 * se],72 "supp_rho_hat_median": float(np.median(rho_hat)),73 "supp_ratio_bias_mc_to_analytical": bias_mc / bias_an,74 "supp_discrepancy_se_units_n_minus_1": (bias_mc - (-(1 + 3 * rho) / (n - 1))) / se,75 "supp_ratio_sd_mc_to_se_asymptotic": sd / float(np.sqrt((1 - rho**2) / n)),76}77with open(os.path.join(ROOT, "out", "results.json"), "w") as f:78 json.dump(res, f, indent=2)79print(json.dumps(res, indent=2))80