import numpy as np, pandas as pd, matplotlib.pyplot as plt, warnings, time
warnings.filterwarnings("ignore")
from sklearn.covariance import LedoitWolf
from scipy.optimize import minimize
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
d=pd.read_csv("stocks_weekly.csv",index_col=0); R=d.values/100.0; N=R.shape[1]; tickers=list(d.columns)
dates=pd.to_datetime(d.index)
print(f"{R.shape[0]} weeks x {N} stocks ({dates.min().date()} to {dates.max().date()})")
print(f"Covariance has N(N+1)/2 = {N*(N+1)//2} free parameters -- vs a 1-year window of only 52 observations.")
print("Tickers:", ", ".join(tickers))

312 weeks x 48 stocks (2019-01-08 to 2024-12-24)
Covariance has N(N+1)/2 = 1176 free parameters -- vs a 1-year window of only 52 observations.
Tickers: AAPL, ABBV, ABT, ADBE, AMD, AMZN, AVGO, AXP, BA, BAC, BLK, C, CAT, COP, COST, CRM, CSCO, CVX, GE, GOOGL, GS, HD, HON, INTC, JNJ, JPM, KO, LLY, LMT, LOW, MCD, META, MRK, MS, MSFT, NKE, NVDA, ORCL, PEP, PFE, PG, TMO, TSLA, UNH, UPS, WFC, WMT, XOM

def minvar(S):
    iv=np.linalg.solve(S, np.ones(N)); return iv/iv.sum()
w_samp=minvar(np.cov(R[:52].T))
fig,ax=plt.subplots(1,2,figsize=(13.5,4.6))
o=np.argsort(w_samp); ax[0].bar(range(N),w_samp[o],color=[RED if v<0 else BLUE for v in w_samp[o]])
ax[0].set_title(f"Sample min-var weights, 1-yr window (gross leverage {np.abs(w_samp).sum():.1f}x)")
ax[0].set_xlabel("stock (sorted)"); ax[0].set_ylabel("portfolio weight"); ax[0].axhline(0,color="k",lw=.6)
# OOS vol vs window length: sample vs Ledoit-Wolf
def roll_vol(win, estim):
    r=[]
    for k in range(win,len(R)-1):
        S=estim(R[k-win:k]); r.append(minvar(S)@R[k])
    return np.std(r)*np.sqrt(52)
# T must EXCEED N: a covariance from T observations has rank <= T-1, so at T <= 48 it is exactly
# singular and the minimum-variance portfolio is undefined. np.linalg.solve does not raise there --
# floating-point noise makes a rank-44 matrix look invertible (condition number ~1e18) and it returns
# garbage. Plotting those points would dramatise the curse with numbers that are not estimates at all.
wins=list(range(N+2,205,15))
vs=[roll_vol(w, lambda X:np.cov(X.T)) for w in wins]
vl=[roll_vol(w, lambda X:LedoitWolf().fit(X).covariance_) for w in wins]
ax[1].plot(wins,np.array(vs)*100,"o-",color=RED,lw=2,label="sample covariance")
ax[1].plot(wins,np.array(vl)*100,"o-",color=BLUE,lw=2,label="Ledoit-Wolf shrinkage")
ax[1].axvline(N,color="k",ls=":",lw=1); ax[1].text(N+2,max(vs)*80,f"N={N} assets",fontsize=8)
ax[1].set_xlabel("estimation window (weeks)"); ax[1].set_ylabel("OOS annualized vol (%)"); ax[1].set_title("The N/T curse: sample risk explodes at short windows"); ax[1].legend()
plt.tight_layout(); plt.show()
c52=np.linalg.cond(np.cov(R[:52].T)); c150=np.linalg.cond(np.cov(R[:150].T))
print(f"At a 1-year window the sample min-var portfolio carries {np.abs(w_samp).sum():.0f}x gross leverage, and its OOS")
print(f"volatility climbs from {vs[-1]*100:.0f}% at the longest window to {vs[0]*100:.0f}% as T approaches N=48, while Ledoit-Wolf stays")
print(f"flat near {np.mean(vl)*100:.0f}% throughout. The mechanism is visible in the conditioning: cond(Sigma) is {c150:.0e} at a")
print(f"150-week window and {c52:.0e} at 52 weeks, and at T<=N it is infinite -- the matrix has no inverse to take.")

At a 1-year window the sample min-var portfolio carries 14x gross leverage, and its OOS
volatility climbs from 12% at the longest window to 61% as T approaches N=48, while Ledoit-Wolf stays
flat near 13% throughout. The mechanism is visible in the conditioning: cond(Sigma) is 9e+02 at a
150-week window and 5e+04 at 52 weeks, and at T<=N it is infinite -- the matrix has no inverse to take.

def longonly_minvar(S):
    r=minimize(lambda w: w@S@w, np.ones(N)/N, method="SLSQP", bounds=[(0,1)]*N,
               constraints=({"type":"eq","fun":lambda w:w.sum()-1},), options={"maxiter":200,"ftol":1e-11})
    return r.x
def ridge_minvar(X): S=np.cov(X.T); return minvar(S+0.1*np.trace(S)/N*np.eye(N))
methods={"sample":lambda X:minvar(np.cov(X.T)),
         "Ledoit-Wolf":lambda X:minvar(LedoitWolf().fit(X).covariance_),
         "ridge (L2)":ridge_minvar,
         "long-only (sparse L1)":lambda X:longonly_minvar(np.cov(X.T)),
         "equal 1/N":lambda X:np.ones(N)/N}
win=52; res={m:{"ret":[],"lev":[],"w":[]} for m in methods}; nz=[]
for k in range(win,len(R)-1):
    X=R[k-win:k]
    for m,f in methods.items():
        w=f(X); res[m]["ret"].append(float(w@R[k])); res[m]["lev"].append(float(np.abs(w).sum())); res[m]["w"].append(w)
    nz.append(int((res["long-only (sparse L1)"]["w"][-1]>1e-4).sum()))
def ann_vol(r): return np.std(r)*np.sqrt(52)*100
def ann_sh(r): return np.mean(r)/np.std(r)*np.sqrt(52)
def turn(ws): W=np.array(ws); return np.mean(np.abs(np.diff(W,axis=0)).sum(1))
tab=pd.DataFrame({m:{"OOS vol %":ann_vol(res[m]["ret"]),"OOS Sharpe":ann_sh(res[m]["ret"]),
                     "gross lev":np.mean(res[m]["lev"]),"turnover":turn(res[m]["w"])} for m in methods}).T
print(tab.round(2).to_string())
print(f"\nlong-only sparse portfolio holds on average {np.mean(nz):.0f} of {N} stocks (the rest get zero weight).")
rngb=np.random.default_rng(0); T=len(res["equal 1/N"]["ret"])
def bootgap(a,b,stat,B=2000):
    A=np.array(res[a]["ret"]); C=np.array(res[b]["ret"])
    g=np.array([stat(A[i])-stat(C[i]) for i in rngb.integers(0,T,(B,T))])
    return g
gs=bootgap("long-only (sparse L1)","equal 1/N",ann_sh)
gv=bootgap("long-only (sparse L1)","equal 1/N",ann_vol)
gl=bootgap("long-only (sparse L1)","Ledoit-Wolf",ann_vol)
d_sh=ann_sh(res["long-only (sparse L1)"]["ret"])-ann_sh(res["equal 1/N"]["ret"])
d_vo=ann_vol(res["long-only (sparse L1)"]["ret"])-ann_vol(res["equal 1/N"]["ret"])
d_lw=ann_vol(res["long-only (sparse L1)"]["ret"])-ann_vol(res["Ledoit-Wolf"]["ret"])
print(f"\nSparse long-only against 1/N, with a paired bootstrap over the {T} weeks rather than an eyeball:")
print(f"   Sharpe     {d_sh:+.3f}  95% CI [{np.percentile(gs,2.5):+.3f}, {np.percentile(gs,97.5):+.3f}]  P(better) = {(gs>0).mean():.2f}")
print(f"   volatility {d_vo:+.2f}pp 95% CI [{np.percentile(gv,2.5):+.2f}, {np.percentile(gv,97.5):+.2f}]")
print(f"   vs Ledoit-Wolf, volatility {d_lw:+.2f}pp 95% CI [{np.percentile(gl,2.5):+.2f}, {np.percentile(gl,97.5):+.2f}]")
print("So: on RISK-ADJUSTED RETURN the sparse portfolio and 1/N are indistinguishable -- the CI straddles zero and the")
print("coin lands heads half the time. On VOLATILITY, which is the quantity these portfolios actually minimise, the")
print("advantage over 1/N is real and the tie with Ledoit-Wolf is genuine. Judge an estimator on its own objective.")

                       OOS vol %  OOS Sharpe  gross lev  turnover
sample                     46.00       -0.07      15.10      7.85
Ledoit-Wolf                15.50        0.36       2.15      0.33
ridge (L2)                 16.26        0.21       2.65      0.48
long-only (sparse L1)      15.00        0.73       1.00      0.20
equal 1/N                  20.49        0.70       1.00      0.00

long-only sparse portfolio holds on average 12 of 48 stocks (the rest get zero weight).

Sparse long-only against 1/N, with a paired bootstrap over the 259 weeks rather than an eyeball:
   Sharpe     +0.028  95% CI [-0.584, +0.606]  P(better) = 0.53
   volatility -5.49pp 95% CI [-8.25, -2.80]
   vs Ledoit-Wolf, volatility -0.51pp 95% CI [-1.51, +0.42]
So: on RISK-ADJUSTED RETURN the sparse portfolio and 1/N are indistinguishable -- the CI straddles zero and the
coin lands heads half the time. On VOLATILITY, which is the quantity these portfolios actually minimise, the
advantage over 1/N is real and the tie with Ledoit-Wolf is genuine. Judge an estimator on its own objective.

fig,ax=plt.subplots(1,3,figsize=(15,4.4)); ms=list(methods)
cols=[RED,BLUE,GREEN,PURP,GREY]
ax[0].bar(ms,[ann_vol(res[m]["ret"]) for m in ms],color=cols); ax[0].set_ylabel("OOS annualized vol (%)"); ax[0].set_title("Out-of-sample volatility (objective)")
ax[1].bar(ms,[np.mean(res[m]["lev"]) for m in ms],color=cols); ax[1].set_ylabel("avg gross leverage"); ax[1].set_title("Leverage (|weights| sum)")
ax[2].bar(ms,[ann_sh(res[m]["ret"]) for m in ms],color=cols); ax[2].set_ylabel("OOS Sharpe"); ax[2].set_title("Risk-adjusted return")
for a in ax: plt.setp(a.get_xticklabels(),rotation=25,ha="right",fontsize=8)
plt.tight_layout(); plt.show()
print("Shrinkage/constraints cut OOS volatility and leverage by an order of magnitude vs sample Markowitz. Ledoit-Wolf and")
print("the sparse long-only portfolio tie for the lowest OOS volatility (~15%); the sparse one uses far fewer holdings/turnover.")

Shrinkage/constraints cut OOS volatility and leverage by an order of magnitude vs sample Markowitz. Ledoit-Wolf and
the sparse long-only portfolio tie for the lowest OOS volatility (~15%); the sparse one uses far fewer holdings/turnover.

wl=longonly_minvar(np.cov(R[-52:].T))
held=np.where(wl>1e-4)[0]; ho=held[np.argsort(-wl[held])]
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
ax[0].bar([tickers[i] for i in ho],[wl[i] for i in ho],color=PURP)
ax[0].set_title(f"Sparse long-only portfolio (latest year): {len(held)}/{N} stocks held"); ax[0].set_ylabel("weight")
plt.setp(ax[0].get_xticklabels(),rotation=90,fontsize=8)
ax[1].hist(nz,bins=range(min(nz),max(nz)+2),color=PURP,edgecolor="white",align="left")
ax[1].axvline(np.mean(nz),color=RED,lw=2,label=f"avg {np.mean(nz):.0f} held"); ax[1].set_xlabel("# stocks held"); ax[1].set_ylabel("weeks"); ax[1].set_title(f"Sparsity over the backtest (of {N} stocks)"); ax[1].legend()
plt.tight_layout(); plt.show()
print(f"The sparse portfolio concentrates in ~{np.mean(nz):.0f} of {N} names and zeroes the rest -- asset selection by L1,")
print("the portfolio counterpart of the factor selection in the previous notebook.")

The sparse portfolio concentrates in ~12 of 48 names and zeroes the rest -- asset selection by L1,
the portfolio counterpart of the factor selection in the previous notebook.

method	OOS vol	leverage	note
sample min-variance	highest (~46%)	~15×	error-maximiser
Ledoit–Wolf shrinkage	≈lowest (~15.5%)	~2×	the asset-risk standard
ridge (L2) covariance	low (~16%)	~2.5×	direct L2 analogue
long-only sparse (L1)	≈lowest (~15%)	1× (no shorts)	selects ~12/48; Sharpe tied with 1/N
equal-weight 1/N	moderate (~20%)	1×	strong Sharpe, no estimation

High-dimensional portfolios — when shrinkage becomes essential¶

Covariance shrinkage & sparse asset selection on 48 stocks, using the Asset-Risk data¶

1. The data and the minimum-variance objective¶

2. The blow-up — sample Markowitz in high dimensions¶

3. The horse race — shrinkage, sparsity, and 1/N¶

4. Reading the race¶

5. The sparse portfolio — asset selection¶

6. Summary¶