""" Q20 — Effective independent names in the 50-ETF book (clean lake). Eigenvalue analysis on the 50-ETF correlation matrix to determine how many effective independent names exist in the book. Output: eigenanalysis.csv + eigenvalue_spectrum.png + stdout summary. """ import pathlib, json import numpy as np import pandas as pd from scipy import linalg LAKE = pathlib.Path("/home/data/lake/market=US/timeframe=1d") OUT = pathlib.Path(__file__).parent 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_returns(start="2026-01-04", end="2026-08-10"): 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) df.columns = [c.lower() for c in df.columns] 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) < 20: continue rets = df[close_col].pct_change().dropna() if len(rets) < 20: continue frames.append(rets.rename(sym)) return pd.DataFrame(frames).T.sort_index() def marchenko_pastur_bound(N, T, q=None): """ Marchenko-Pastur upper bound for eigenvalues of a random correlation matrix. q = T/N ratio. Eigenvalues above this bound are 'signal'. """ if q is None: q = T / N sigma2 = 1.0 # correlation matrix has unit diagonal lambda_plus = sigma2 * (1 + 1/np.sqrt(q))**2 return lambda_plus def participation_ratio(eigenvalues): """Participation ratio: (sum(lambda))^2 / sum(lambda^2). Equals N for identity.""" lam = eigenvalues[eigenvalues > 0] return (np.sum(lam))**2 / np.sum(lam**2) def main(): print("Loading 50-ETF daily returns (test window: 2026-01-04 to 2026-08-10)...") rets = load_etf_returns() N = rets.shape[0] # symbols (rows) T = rets.shape[1] # trading days (columns) print(f"Loaded {N} ETFs, {T} trading days") print(f"Note: N={N} symbols (rows), T={T} days (columns) in return matrix") # Drop any ETFs with too many NaNs rets = rets.dropna(axis=0, thresh=int(T * 0.8)) N = rets.shape[0] rets = rets.fillna(0) print(f"After dropping high-NaN ETFs: {N} symbols") # Correlation matrix corr = rets.T.corr() print(f"Correlation matrix: {corr.shape}") # Eigendecomposition eigvals_raw = linalg.eigvalsh(corr.values) eigvals = np.sort(eigvals_raw)[::-1] # descending # Marchenko-Pastur bound q_ratio = T / N mp_bound = marchenko_pastur_bound(N, T, q_ratio) n_signal = int(np.sum(eigvals > mp_bound)) print(f"\n=== Eigenvalue Analysis ===") print(f" N (ETFs): {N}") print(f" T (days): {T}") print(f" q = T/N: {q_ratio:.2f}") print(f" Marchenko-Pastur upper bound: {mp_bound:.4f}") print(f" Eigenvalues above MP bound (signal): {n_signal}") print(f"\n Top 10 eigenvalues:") for i, ev in enumerate(eigvals[:10]): pct = ev / eigvals.sum() * 100 marker = " * SIGNAL" if ev > mp_bound else "" print(f" λ_{i+1:2d} = {ev:8.4f} ({pct:5.1f}% var){marker}") # Cumulative variance share cumvar = np.cumsum(eigvals) / eigvals.sum() print(f"\n Cumulative variance explained by top-k components:") for k in [1, 2, 3, 4, 5, 10, 15, 20]: if k <= len(cumvar): print(f" Top {k:2d}: {cumvar[k-1]*100:5.1f}%") # Effective rank measures pr = participation_ratio(eigvals) # 80% variance count var_80 = int(np.searchsorted(cumvar, 0.80) + 1) # 90% variance count var_90 = int(np.searchsorted(cumvar, 0.90) + 1) print(f"\n Participation ratio (effective rank): {pr:.2f}") print(f" Eigenvalues needed for 80% variance: {var_80}") print(f" Eigenvalues needed for 90% variance: {var_90}") # --- Save --- eigen_df = pd.DataFrame({ "rank": range(1, len(eigvals) + 1), "eigenvalue": eigvals, "pct_variance": eigvals / eigvals.sum() * 100, "cumulative_pct": cumvar * 100, "above_mp_bound": eigvals > mp_bound, }) eigen_df.to_csv(OUT / "eigenanalysis.csv", index=False) summary = { "N_symbols": N, "T_days": T, "q_ratio": round(q_ratio, 2), "mp_bound": round(float(mp_bound), 4), "n_signal_eigenvalues": n_signal, "top_eigenvalues": [round(float(ev), 4) for ev in eigvals[:10]], "top_pct_variance": [round(float(ev / eigvals.sum() * 100), 1) for ev in eigvals[:10]], "participation_ratio": round(float(pr), 2), "eigenvalues_for_80pct_var": var_80, "eigenvalues_for_90pct_var": var_90, "cumulative_var_top4": round(float(cumvar[3] * 100), 1) if len(cumvar) > 3 else None, "cumulative_var_top10": round(float(cumvar[9] * 100), 1) if len(cumvar) > 9 else None, } with open(OUT / "eigen_summary.json", "w") as f: json.dump(summary, f, indent=2) print(f"\nSaved: {OUT / 'eigenanalysis.csv'}") print(f"Saved: {OUT / 'eigen_summary.json'}") # --- Verdict --- print(f"\n=== VERDICT ===") if var_80 <= 5: print(f" CONFIRMED: top-{var_80} components explain 80%+ of variance.") print(f" The 50-ETF book has ≈{var_80} effective independent names.") print(f" This explains why topk 10→20 adds no breadth (EVIDENCE#024).") elif var_80 <= 10: print(f" PARTIAL: top-{var_80} for 80% variance — moderate concentration.") print(f" Participation ratio = {pr:.1f}, suggesting ~{pr:.0f} effective names.") else: print(f" REFUTED: need {var_80} components for 80% variance — book is well-diversified.") print(f" The 'only ~4 effective names' claim is overstated.") if __name__ == "__main__": main()