166 lines
6.1 KiB
Python
166 lines
6.1 KiB
Python
"""
|
|
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()
|