1"""F5-conf-02: finite-sample bias of OLS AR(1) coefficient (with intercept).2Reads only data/design.csv; writes out/results.json deterministically.3Run from workspace root: ~/wp29/venv/bin/python out/analysis.py"""4import json, math, os5import numpy as np6import pandas as pd78ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))9design = pd.read_csv(os.path.join(ROOT, "data", "design.csv"))10p = dict(zip(design["parameter"], design["value"]))11data_issues = []1213n = 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) < 1 and sigma > 0 and n >= 3 and R >= 21819# Baseline / analytical quantities20bias_analytical = -(1 + 3 * rho) / n21se_asymptotic = math.sqrt((1 - rho ** 2) / n)2223# Monte Carlo (vectorised across replications); c = 0 in DGP24rng = np.random.default_rng(seed)25y = np.empty((R, n))26y[:, 0] = rng.normal(0.0, sigma / math.sqrt(1 - rho ** 2), size=R)27eps = rng.normal(0.0, sigma, size=(R, n - 1))28for t in range(1, n):29 y[:, t] = rho * y[:, t - 1] + eps[:, t - 1]3031x = y[:, :-1]; z = y[:, 1:] # regress y_t on (1, y_{t-1}), t = 2..n32xd = x - x.mean(axis=1, keepdims=True)33zd = z - z.mean(axis=1, keepdims=True)34rho_hat = (xd * zd).sum(axis=1) / (xd * xd).sum(axis=1)3536finite = np.isfinite(rho_hat)37if not finite.all():38 data_issues.append({"issue": f"{(~finite).sum()} non-finite rho_hat", "handling": "excluded"})39rh = rho_hat[finite]40R_used = int(rh.size)4142rho_hat_mean = float(rh.mean())43rho_hat_sd = float(rh.std(ddof=1))44bias_mc = rho_hat_mean - rho45bias_mc_se = rho_hat_sd / math.sqrt(R_used)46discrepancy = (bias_mc - bias_analytical) / bias_mc_se4748# Sensitivity (not a required key): formula with effective sample size n-149bias_analytical_nm1 = -(1 + 3 * rho) / (n - 1)5051results = {52 "task_id": "F5-conf-02",53 "bias_analytical": bias_analytical,54 "bias_mc": bias_mc,55 "bias_mc_se": bias_mc_se,56 "rho_hat_mean": rho_hat_mean,57 "rho_hat_sd_mc": rho_hat_sd,58 "discrepancy_se_units": discrepancy,59 "se_asymptotic": se_asymptotic,60 "n_replications": R,61 "feasibility": "complete",62 "feasibility_reason": None,63 "data_issues": data_issues,64 "n_input_rows": {"design.csv": int(len(design))},65 "n_analysed": R_used,66 "random_seed": seed,67 # supplementary68 "baseline_naive_bias": 0.0,69 "bias_analytical_effective_n_minus_1": bias_analytical_nm1,70 "discrepancy_se_units_n_minus_1": (bias_mc - bias_analytical_nm1) / bias_mc_se,71 "rho_hat_median": float(np.median(rh)),72 "rho_hat_sd_ratio_to_asymptotic": rho_hat_sd / se_asymptotic,73 "relative_error_analytical_vs_mc": (bias_analytical - bias_mc) / bias_mc,74 "design": {"n": n, "rho": rho, "sigma": sigma, "intercept": True,75 "initial_condition": "stationary", "c": 0.0},76}77os.makedirs(os.path.join(ROOT, "out"), exist_ok=True)78with open(os.path.join(ROOT, "out", "results.json"), "w") as f:79 json.dump(results, f, indent=2)80for k in ["bias_analytical", "bias_mc", "bias_mc_se", "discrepancy_se_units",81 "bias_analytical_effective_n_minus_1", "discrepancy_se_units_n_minus_1",82 "rho_hat_mean", "rho_hat_median", "rho_hat_sd_mc", "se_asymptotic",83 "relative_error_analytical_vs_mc"]:84 print(f"{k}: {results[k]:.6f}")85