218 lines
8.6 KiB
Python
218 lines
8.6 KiB
Python
"""
|
||
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()
|