import numpy as np, pandas as pd, matplotlib.pyplot as plt, warnings
warnings.filterwarnings("ignore")
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
p=pd.read_csv("factor_zoo.csv",index_col="ym"); fac=[c for c in p.columns if c!="RF"]
R=p[fac]; N=len(fac)
dates=pd.to_datetime(p.index.astype(str),format="%Y%m")
prem=R.mean()*12; sharpe1=R.mean()/R.std()*np.sqrt(12)
print(f"{N} factors, {len(R)} months, {p.index.min()//100}-{p.index.max()//100}")
fig,ax=plt.subplots(figsize=(11,4.2))
o=prem.sort_values(); ax.bar(range(N),o.values,color=[GREEN if v>0 else RED for v in o.values])
ax.set_xticks(range(N)); ax.set_xticklabels(o.index,rotation=90,fontsize=8); ax.axhline(0,color="k",lw=.7)
ax.set_ylabel("annualized mean premium (%)"); ax.set_title("The factor zoo — 20 long-short factors (1963-2026)")
plt.tight_layout(); plt.show()
print("Positive-premium factors pay to go long the high leg (Mom, HML, RMW...); negative ones pay to SHORT it")
print("(high investment, high issuance, high volatility underperform). The SDF can weight either direction.")

20 factors, 755 months, 1963-2026

Positive-premium factors pay to go long the high leg (Mom, HML, RMW...); negative ones pay to SHORT it
(high investment, high issuance, high volatility underperform). The SDF can weight either direction.

fig,ax=plt.subplots(1,2,figsize=(13,5))
C=R.corr().values
im=ax[0].imshow(C,cmap="RdBu_r",vmin=-1,vmax=1); ax[0].set_xticks(range(N)); ax[0].set_xticklabels(fac,rotation=90,fontsize=7)
ax[0].set_yticks(range(N)); ax[0].set_yticklabels(fac,fontsize=7); ax[0].set_title("Factor correlation matrix — themes clump")
plt.colorbar(im,ax=ax[0],fraction=0.046)
ev=np.linalg.eigvalsh(np.corrcoef((R-R.mean()).T))[::-1]; cum=np.cumsum(ev)/ev.sum()
ax[1].bar(range(1,N+1),ev,color=GREY); ax[1].set_xlabel("principal component"); ax[1].set_ylabel("eigenvalue",color=GREY)
a2=ax[1].twinx(); a2.plot(range(1,N+1),cum*100,"o-",color=BLUE,lw=2); a2.set_ylabel("cumulative variance (%)",color=BLUE)
a2.axhline(90,color=RED,ls=":",lw=1); k90=np.argmax(cum>=0.9)+1; a2.text(N*0.5,80,f"{k90} PCs explain 90%",color=RED,fontsize=9)
ax[1].set_title("PCA scree — few independent dimensions")
plt.tight_layout(); plt.show()
rp=np.random.default_rng(0)
def _k90(M):
    e=np.linalg.eigvalsh(np.corrcoef(M.T))[::-1]; return int(np.argmax(np.cumsum(e)/e.sum()>=0.9))+1
knull=np.mean([_k90(rp.standard_normal((len(R),N))) for _ in range(200)])
print(f"{k90} principal components are needed to reach 90% of the covariation among {N} factors. The benchmark that")
print(f"makes that number mean something: {N} INDEPENDENT series of the same length need {knull:.0f}. So the zoo does carry")
print(f"real redundancy -- but not the one-or-two-themes story sometimes told: PC1 alone is {100*cum[0]:.0f}%, and it still")
print(f"takes {int(np.argmax(cum>=0.75))+1} components to reach 75%. Enough duplication for a selector to exploit, not enough to call it a facade.")

10 principal components are needed to reach 90% of the covariation among 20 factors. The benchmark that
makes that number mean something: 20 INDEPENDENT series of the same length need 18. So the zoo does carry
real redundancy -- but not the one-or-two-themes story sometimes told: PC1 alone is 32%, and it still
takes 5 components to reach 75%. Enough duplication for a selector to exploit, not enough to call it a facade.

from sklearn.linear_model import LinearRegression, RidgeCV, LassoCV, ElasticNetCV, lasso_path
tr=p.index<=200012; te=p.index>200012
Rtr=R.values[tr]; Rte=R.values[te]; sd=Rtr.std(0); Ztr=Rtr/sd; Zte=Rte/sd; y=np.ones(len(Ztr))
def sharpe(x): return np.sqrt(12)*x.mean()/x.std()
bols=LinearRegression(fit_intercept=False).fit(Ztr,y).coef_
is_s=sharpe(Ztr@bols); oos_s=sharpe(Zte@bols)
print(f"Kitchen-sink tangency SDF ({N} factors):")
print(f"   IN-SAMPLE  Sharpe (1963-2000): {is_s:5.2f}   <- the mirage")
print(f"   OUT-OF-SAMPLE Sharpe (2001-26): {oos_s:5.2f}   <- the reality")
print(f"   best SINGLE factor OOS Sharpe : {max(sharpe(Zte[:,j]) for j in range(N)):5.2f}")
fig,ax=plt.subplots(figsize=(6.2,4.2))
ax.bar(["in-sample\n(1963-2000)","out-of-sample\n(2001-2026)"],[is_s,oos_s],color=[GREY,RED])
for i,v in enumerate([is_s,oos_s]): ax.text(i,v+0.03,f"{v:.2f}",ha="center")
ax.set_ylabel("annualized Sharpe"); ax.set_title("Kitchen-sink SDF: in-sample mirage vs OOS reality")
plt.tight_layout(); plt.show()
best_j=int(np.argmax([sharpe(Zte[:,j]) for j in range(N)]))
print("A Sharpe of ~2.4 in-sample is not skill -- it is Sigma^{-1}mu fitting noise across 20 correlated factors.")
print(f"Worth pausing on the third number: the 20-factor tangency portfolio ({oos_s:.2f}) is BEATEN out of sample by simply")
print(f"holding {fac[best_j]} on its own ({sharpe(Zte[:,best_j]):.2f}). Sixty years of data, twenty published factors, and the optimiser")
print("cannot outperform one of its own inputs -- that is how much damage estimation error does to Sigma^{-1}mu.")

Kitchen-sink tangency SDF (20 factors):
   IN-SAMPLE  Sharpe (1963-2000):  2.41   <- the mirage
   OUT-OF-SAMPLE Sharpe (2001-26):  0.40   <- the reality
   best SINGLE factor OOS Sharpe :  0.55

A Sharpe of ~2.4 in-sample is not skill -- it is Sigma^{-1}mu fitting noise across 20 correlated factors.
Worth pausing on the third number: the 20-factor tangency portfolio (0.40) is BEATEN out of sample by simply
holding Mkt-RF on its own (0.55). Sixty years of data, twenty published factors, and the optimiser
cannot outperform one of its own inputs -- that is how much damage estimation error does to Sigma^{-1}mu.

alphas,coefs,_=lasso_path(Ztr,y,n_alphas=60,eps=1e-3)
nfac=(np.abs(coefs)>1e-8).sum(0); is_path=[sharpe(Ztr@coefs[:,k]) for k in range(coefs.shape[1])]
oos_path=[sharpe(Zte@coefs[:,k]) if np.any(coefs[:,k]) else 0 for k in range(coefs.shape[1])]
lc=LassoCV(n_alphas=100,cv=5,fit_intercept=False,max_iter=100000).fit(Ztr,y)
b_lasso=lc.coef_; k_l=int((np.abs(b_lasso)>1e-8).sum())
fig,ax=plt.subplots(1,2,figsize=(13.5,4.8))
ax[0].plot(nfac,is_path,"o-",color=GREY,lw=2,label="in-sample Sharpe")
ax[0].plot(nfac,oos_path,"o-",color=BLUE,lw=2,label="out-of-sample Sharpe")
ax[0].axvline(k_l,color=GREEN,ls="--",label=f"Lasso CV pick: {k_l} factors")
ax[0].set_xlabel("number of factors in the SDF"); ax[0].set_ylabel("annualized Sharpe"); ax[0].set_title("Overfitting scissors: fit rises, OOS plateaus"); ax[0].legend(fontsize=8)
srt=np.argsort(np.abs(b_lasso)); srt=srt[np.abs(b_lasso[srt])>1e-8]
ax[1].barh([fac[i] for i in srt],[b_lasso[i] for i in srt],color=[GREEN if b_lasso[i]>0 else RED for i in srt])
ax[1].set_xlabel("SDF weight (standardized factors)"); ax[1].set_title(f"The sparse SDF: {k_l} factors Lasso keeps")
plt.tight_layout(); plt.show()
oa=np.array(oos_path,float); pk=int(np.nanargmax(oa))
VAL=["HML","Value_BEME","Earn_EP","Cash_CFP","Div_DP"]; LOW=["Var","ResVar","Beta"]
kept=[fac[i] for i in srt[::-1]]
print(f"Lasso-CV SDF: {k_l} of {N} factors, OOS Sharpe {sharpe(Zte@b_lasso):.2f} against the kitchen sink's {oos_s:.2f} -- it does not")
print(f"merely match the 20-factor fit, it edges it, and the path peaks at {oa[pk]:.2f} with {nfac[pk]} factors. Cross-validation lands")
print("essentially on that peak without ever seeing the test period, which is the result worth having.")
print(f"Kept: {', '.join(kept)}.")
print(f"The pruning is exactly where the redundancy is: {sum(f in VAL for f in kept)} of the {len(VAL)} value proxies survives "
      f"({', '.join(f for f in kept if f in VAL)}) and {sum(f in LOW for f in kept)} of the {len(LOW)} low-risk proxies "
      f"({', '.join(f for f in kept if f in LOW)}).")

Lasso-CV SDF: 10 of 20 factors, OOS Sharpe 0.44 against the kitchen sink's 0.40 -- it does not
merely match the 20-factor fit, it edges it, and the path peaks at 0.45 with 12 factors. Cross-validation lands
essentially on that peak without ever seeing the test period, which is the result worth having.
Kept: Mom, ResVar, ST_Rev, SMB, Mkt-RF, LT_Rev, RMW, CMA, HML, Accr_AC.
The pruning is exactly where the redundancy is: 1 of the 5 value proxies survives (HML) and 1 of the 3 low-risk proxies (ResVar).

from sklearn.linear_model import RidgeCV
def norm(b): s=np.sum(np.abs(b)); return b/s if s>0 else b
w_ks=norm(bols); w_l=norm(b_lasso); w_eq=norm(np.sign(Rtr.mean(0)))
ret={"kitchen-sink OLS":Zte@w_ks,"Lasso sparse SDF":Zte@w_l,"equal-weight (signs)":Zte@w_eq}
dte=dates[te]
fig,ax=plt.subplots(figsize=(9,4.4))
for (nm,r_),c in zip(ret.items(),[GREY,BLUE,ORANGE]):
    ax.plot(dte,np.cumsum(r_),color=c,lw=2,label=f"{nm}  (OOS Sharpe {sharpe(r_):.2f})")
ax.set_ylabel("cumulative OOS return (%)"); ax.set_title("Out-of-sample SDF performance (2001-2026)"); ax.legend(fontsize=9); ax.grid(alpha=.3)
plt.tight_layout(); plt.show()
print("Verdict: the sparse Lasso SDF matches the 20-factor kitchen sink out of sample with ~half the factors -- parsimony")
print("at no cost, and a far more interpretable model. It does NOT dramatically beat OLS on Sharpe here: with 60 years of")
print("20 well-diversified factors the sample tangency is already decent. Shrinkage's OOS edge grows with the DIMENSION of")
print("the problem -- Kozak-Nagel-Santosh need hundreds of characteristic portfolios for it to dominate.")

Verdict: the sparse Lasso SDF matches the 20-factor kitchen sink out of sample with ~half the factors -- parsimony
at no cost, and a far more interpretable model. It does NOT dramatically beat OLS on Sharpe here: with 60 years of
20 well-diversified factors the sample tangency is already decent. Shrinkage's OOS edge grows with the DIMENSION of
the problem -- Kozak-Nagel-Santosh need hundreds of characteristic portfolios for it to dominate.

Finance application — shrinking the factor zoo¶

Lasso / Elastic-Net factor selection in the cross-section of returns¶

1. The factor zoo — 20 long-short factors¶

What a "factor" actually is¶

2. The zoo is redundant¶

3. The in-sample mirage¶

4. Lasso selection — taming the zoo¶

5. Out-of-sample, and the honest verdict¶

6. Summary¶