""" Q19 — Variance-ratio study on the 50-ETF panel (clean lake). Tests whether assets are submartingales long-horizon / mean-reverting short-horizon (VR < 1 at 5–20d). Uses the Lo–MacKinlay heteroskedasticity- robust VR statistic. Output: VR_stats.csv + stdout summary. """ import pathlib, json, sys import numpy as np import pandas as pd from scipy import stats LAKE = pathlib.Path("/home/data/lake/market=US/timeframe=1d") OUT = pathlib.Path(__file__).parent # --- 50-ETF panel (all non-single-stock names in the lake) --- SINGLE_STOCKS = { "AAPL","MSFT","NVDA","AMZN","GOOGL","META","TSLA","AVGO","AMD", "JPM","UNH","PG","JNJ","MA","V","WMT","DIS","HD","KO","PEP", "BAC","XOM","MCD","ABBV","COST","CRM","NFLX","ORCL","IBM","T", } def load_etf_bars(start="2015-01-01", end="2026-08-19"): frames = [] for f in sorted(LAKE.glob("symbol=*.parquet")): sym = f.stem.replace("symbol=", "") if sym in SINGLE_STOCKS: continue df = pd.read_parquet(f) if len(df) < 100: continue df.columns = [c.lower() for c in df.columns] # Lake uses 'c' for close, 'date' column for date close_col = "c" if "c" in df.columns else "close" if close_col not in df.columns: continue if "date" in df.columns: df = df.set_index("date") elif "datetime" in df.columns: df = df.set_index("datetime") df.index = pd.to_datetime(df.index) df = df.loc[start:end] if len(df) < 200: continue frames.append(df[close_col].rename(sym)) return pd.DataFrame(frames).T.sort_index() def variance_ratio(series, q): """ Lo-MacKinlay variance ratio with heteroskedasticity-robust z-stat. VR(q) = Var(q-period returns) / (q * Var(1-period returns)) H0: VR = 1 (random walk). VR < 1 => mean reversion; VR > 1 => momentum / trending. """ y = series.dropna().values n = len(y) if n < q + 10: return np.nan, np.nan, np.nan rets = np.diff(np.log(y)) n_ret = len(rets) mu = np.mean(rets) # 1-period variance (with heteroskedasticity correction) m2 = np.sum((rets - mu) ** 2) / (n_ret - 1) # q-period returns rq = np.array([np.sum(rets[i:i+q]) for i in range(n_ret - q + 1)]) vq = np.var(rq, ddof=1) vr = vq / (q * m2) if m2 > 0 else np.nan # Robust z-stat (heteroskedasticity-robust, Lo-MacKinlay 1988 Eq. 18) # Under H0: VR=1, z ~ N(0,1) T = n_ret # Sum of autocovariances for q-period returns mu_q = np.mean(rq) # Omega_1 (heteroskedasticity-robust variance of VR estimate) # Simplified: use the asymptotic variance under heteroskedasticity delta = np.zeros(q) for j in range(1, q): rho_j = np.corrcoef(rets[j:], rets[:-j])[0, 1] if len(rets) > j + 1 else 0 delta[j] = 2 * (1 - j/q) * rho_j omega2 = np.sum(delta) # z-stat se_vr = np.sqrt(max((2 * (2*q - 1) * (q-1)) / (3 * q * T) * (1 + omega2), 1e-15)) z = (vr - 1) / se_vr if se_vr > 0 else 0 pval = 2 * (1 - stats.norm.cdf(abs(z))) return vr, z, pval def main(): print("Loading 50-ETF daily bars from lake...") prices = load_etf_bars() print(f"Loaded {prices.shape[1]} symbols, {prices.shape[0]} trading days ({prices.index[0].date()} to {prices.index[-1].date()})") horizons = [5, 10, 20] results = [] for sym in prices.columns: s = prices[sym].dropna() if len(s) < 500: continue row = {"symbol": sym, "n_days": len(s)} for q in horizons: vr, z, p = variance_ratio(s, q) row[f"VR_{q}d"] = round(vr, 4) row[f"z_{q}d"] = round(z, 2) row[f"p_{q}d"] = round(p, 4) results.append(row) df = pd.DataFrame(results) # --- Summary --- print("\n=== Variance Ratio Summary (50-ETF Panel, 2015-01-01 to 2026-08-19) ===") for q in horizons: vr_col = f"VR_{q}d" valid = df[vr_col].dropna() frac_lt1 = (valid < 1).mean() frac_sig_revert = ((valid < 1) & (df[f"z_{q}d"].abs() > 2)).mean() frac_sig_momentum = ((valid > 1) & (df[f"z_{q}d"].abs() > 2)).mean() print(f"\n Horizon {q}d:") print(f" Mean VR: {valid.mean():.4f}, Median VR: {valid.median():.4f}") print(f" Std VR: {valid.std():.4f}") print(f" Fraction VR < 1: {frac_lt1:.1%} ({(valid < 1).sum()}/{len(valid)})") print(f" Fraction VR < 1 & |z|>2 (mean-revert): {frac_sig_revert:.1%}") print(f" Fraction VR > 1 & |z|>2 (momentum): {frac_sig_momentum:.1%}") print(f" Min VR: {valid.min():.4f}, Max VR: {valid.max():.4f}") # --- Cross-check: pooled trend-slope beta --- print("\n=== Cross-check: sp_trend_slope_5 regression ===") # Compute log-price momentum slope for each symbol betas = [] for sym in prices.columns: s = prices[sym].dropna() if len(s) < 100: continue logp = np.log(s.values) # 5-day rolling slope (regress logp on [0,1,2,3,4] for each window) slopes = [] for i in range(len(logp) - 4): y_win = logp[i:i+5] x_win = np.arange(5) # OLS slope slope = (5 * np.sum(x_win * y_win) - np.sum(x_win) * np.sum(y_win)) / (5 * np.sum(x_win**2) - np.sum(x_win)**2) slopes.append(slope) # Future 5-day return rets_5d = np.array([np.log(s.values[i+5] / s.values[i]) for i in range(len(s) - 5)]) slopes_arr = np.array(slopes[:len(rets_5d)]) if len(slopes_arr) < 50: continue # Regression: future 5d return ~ beta * trend_slope_5 valid_mask = np.isfinite(slopes_arr) & np.isfinite(rets_5d) if valid_mask.sum() < 50: continue slope_valid = slopes_arr[valid_mask] ret_valid = rets_5d[valid_mask] # OLS X = np.column_stack([np.ones(len(slope_valid)), slope_valid]) beta_hat = np.linalg.lstsq(X, ret_valid, rcond=None)[0] betas.append({"symbol": sym, "beta": beta_hat[1], "n": valid_mask.sum()}) beta_df = pd.DataFrame(betas) if len(beta_df) > 0: pooled_beta = beta_df["beta"].mean() pooled_se = beta_df["beta"].std() / np.sqrt(len(beta_df)) t_stat = pooled_beta / pooled_se if pooled_se > 0 else 0 print(f" Panel ({len(beta_df)} symbols): mean slope-beta = {pooled_beta:.4f}, SE = {pooled_se:.4f}, t = {t_stat:.2f}") print(f" Beta range: [{beta_df['beta'].min():.4f}, {beta_df['beta'].max():.4f}]") n_negative = (beta_df["beta"] < 0).sum() print(f" Symbols with negative beta (mean-revert): {n_negative}/{len(beta_df)} ({n_negative/len(beta_df):.1%})") # --- Save --- df.to_csv(OUT / "VR_stats.csv", index=False) summary = { "panel_size": len(df), "date_range": f"{prices.index[0].date()} to {prices.index[-1].date()}", "n_trading_days": len(prices), "horizons": {}, } for q in horizons: valid = df[f"VR_{q}d"].dropna() summary["horizons"][f"{q}d"] = { "mean_vr": round(float(valid.mean()), 4), "median_vr": round(float(valid.median()), 4), "frac_lt1": round(float((valid < 1).mean()), 3), "frac_sig_revert_z2": round(float(((valid < 1) & (df[f"z_{q}d"].abs() > 2)).mean()), 3), "frac_sig_momentum_z2": round(float(((valid > 1) & (df[f"z_{q}d"].abs() > 2)).mean()), 3), } if len(beta_df) > 0: summary["trend_slope_5_beta"] = { "mean": round(float(pooled_beta), 4), "se": round(float(pooled_se), 4), "t_stat": round(float(t_stat), 2), "n_negative": int(n_negative), "n_total": len(beta_df), } with open(OUT / "VR_summary.json", "w") as f: json.dump(summary, f, indent=2) print(f"\nSaved: {OUT / 'VR_stats.csv'}") print(f"Saved: {OUT / 'VR_summary.json'}") # --- Verdict --- print("\n=== VERDICT ===") vr5 = summary["horizons"]["5d"] vr10 = summary["horizons"]["10d"] vr20 = summary["horizons"]["20d"] any_revert = any(h["frac_sig_revert_z2"] > 0.1 for h in [vr5, vr10, vr20]) all_lt1_median = all(h["median_vr"] < 1 for h in [vr5, vr10, vr20]) if all_lt1_median and any_revert: print(" SUPPORTS mean-reversion hypothesis: median VR < 1 at all horizons,") print(" material fraction with significant mean-reversion (|z| > 2).") elif all_lt1_median: print(" PARTIAL: median VR < 1 at all horizons, but few significant z-stats.") else: print(" REFUTES strict mean-reversion: median VR >= 1 at some horizons.") print(" See VR_stats.csv for per-symbol detail.") if __name__ == "__main__": main()