import numpy as np, pandas as pd, matplotlib.pyplot as plt, warnings
warnings.filterwarnings("ignore")
from matplotlib.patches import Ellipse
from sklearn.preprocessing import StandardScaler
from cluster import kmeans, GMM
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
d=pd.read_csv("spx_rv_ret.csv"); dates=pd.to_datetime(d["date"]); ret=d["ret"].values; lvol=np.log(np.sqrt(d["rv"].values)*100)
X=np.column_stack([ret,lvol]); scaler=StandardScaler().fit(X); Z=scaler.transform(X)
print(f"{len(X)} trading days, 2 features (daily return %, log realized vol); no labels -- clustering must find the structure")
fig,ax=plt.subplots(figsize=(6.5,5))
ax.scatter(ret,lvol,s=6,alpha=.25,color=GREY); ax.set_xlabel("daily return (%)"); ax.set_ylabel("log realized volatility"); ax.set_title("The market in the return-volatility plane")
plt.tight_layout(); plt.show()
print("A calm low-volatility core plus a high-volatility wing holding both crashes and sharp rebounds -- structure a")
print("clustering algorithm should recover as distinct market regimes.")

3459 trading days, 2 features (daily return %, log realized vol); no labels -- clustering must find the structure

A calm low-volatility core plus a high-volatility wing holding both crashes and sharp rebounds -- structure a
clustering algorithm should recover as distinct market regimes.

from sklearn.cluster import KMeans
lab,C,inertia=kmeans(Z,3,seed=0)
sk=KMeans(3,n_init=10,random_state=0).fit(Z)
print(f"k-means inertia (within-cluster SSE): from-scratch {inertia:.1f}  vs  scikit-learn {sk.inertia_:.1f}  -> match")
Cx=scaler.inverse_transform(C)                                 # centroids back to original units
fig,ax=plt.subplots(figsize=(6.5,5))
for j in range(3): ax.scatter(ret[lab==j],lvol[lab==j],s=6,alpha=.35,label=f"cluster {j}")
ax.scatter(Cx[:,0],Cx[:,1],c="black",marker="X",s=140,label="centroids")
ax.set_xlabel("daily return (%)"); ax.set_ylabel("log realized volatility"); ax.set_title("k-means (k=3): hard, spherical clusters"); ax.legend(fontsize=8)
plt.tight_layout(); plt.show()
print("Inertia matches the reference solver. But k-means draws straight boundaries and forces hard membership -- it cannot")
print("model the elliptical, overlapping shape of volatile-up vs volatile-down days, which is where the GMM comes in.")

k-means inertia (within-cluster SSE): from-scratch 2982.9  vs  scikit-learn 2983.0  -> match

Inertia matches the reference solver. But k-means draws straight boundaries and forces hard membership -- it cannot
model the elliptical, overlapping shape of volatile-up vs volatile-down days, which is where the GMM comes in.

from sklearn.mixture import GaussianMixture
g=GMM(3,seed=0).fit(Z); sg=GaussianMixture(3,covariance_type="full",random_state=0,n_init=5).fit(Z)
print(f"GMM log-likelihood: from-scratch {g.loglik_[-1]:.1f}  vs  scikit-learn {sg.score(Z)*len(Z):.1f}")
print(f"GMM BIC:            from-scratch {g.bic(Z):.1f}  vs  scikit-learn {sg.bic(Z):.1f}")
R=g.predict_proba(Z); lab_g=R.argmax(1); conf=R.max(1)
def ell(ax,mu,Sig,color):
    v,w=np.linalg.eigh(Sig); ang=np.degrees(np.arctan2(w[1,0],w[0,0]))
    for s in (1,2): ax.add_patch(Ellipse(mu,2*s*np.sqrt(v[0]),2*s*np.sqrt(v[1]),angle=ang,fill=False,color=color,lw=1.6,alpha=.8))
fig,ax=plt.subplots(1,2,figsize=(13,4.8)); cols=[BLUE,RED,GREEN]
for j in range(3):
    m=lab_g==j; ax[0].scatter(ret[m],lvol[m],s=6,alpha=.35,color=cols[j])
    mu_o=scaler.inverse_transform(g.mu[j:j+1])[0]; Sig_o=np.diag(scaler.scale_)@g.Sig[j]@np.diag(scaler.scale_)
    ell(ax[0],mu_o,Sig_o,cols[j])
ax[0].set_xlabel("daily return (%)"); ax[0].set_ylabel("log realized volatility"); ax[0].set_title("GMM: elliptical clusters (1 & 2 std ellipses)")
sc=ax[1].scatter(ret,lvol,s=8,c=conf,cmap="viridis",vmin=.5,vmax=1); plt.colorbar(sc,ax=ax[1],label="assignment confidence (max responsibility)")
ax[1].set_xlabel("daily return (%)"); ax[1].set_ylabel("log realized volatility"); ax[1].set_title("Soft assignments: confident in cores, uncertain at borders")
plt.tight_layout(); plt.show()
from sklearn.metrics import adjusted_rand_score
print(f"\nMatching the fit is not the same as matching the PARTITION. Agreement between the two 3-cluster labelings")
print(f"(adjusted Rand index, 1 = identical): {adjusted_rand_score(lab_g,sg.predict(Z)):.3f} -- the same taxonomy, but the boundaries between the two")
print(f"turbulent clusters sit in different places. EM is a local optimizer and the likelihood surface has several modes;")
print(f"here the from-scratch run reaches the higher likelihood of the two, so it is not the weaker solution.")
print("\nThe GMM's ellipses tilt to capture the volatile-up and volatile-down wings, and its SOFT assignments admit")
print("uncertainty on the boundary days -- exactly the days a hard label would misrepresent.")

GMM log-likelihood: from-scratch -8728.2  vs  scikit-learn -8732.1
GMM BIC:            from-scratch 17594.9  vs  scikit-learn 17602.7

Matching the fit is not the same as matching the PARTITION. Agreement between the two 3-cluster labelings
(adjusted Rand index, 1 = identical): 0.896 -- the same taxonomy, but the boundaries between the two
turbulent clusters sit in different places. EM is a local optimizer and the likelihood surface has several modes;
here the from-scratch run reaches the higher likelihood of the two, so it is not the weaker solution.

The GMM's ellipses tilt to capture the volatile-up and volatile-down wings, and its SOFT assignments admit
uncertainty on the boundary days -- exactly the days a hard label would misrepresent.

KS=list(range(1,13))                                           # run the grid well past the candidate answers
inertias=[kmeans(Z,k,seed=0)[2] for k in KS]; bics=[GMM(k,seed=0).fit(Z).bic(Z) for k in KS]
kbic=KS[int(np.argmin(bics))]

# Stability: refit on two 80% subsamples and compare the labels they give the days both saw.
SK=list(range(2,9)); rng=np.random.default_rng(0); stab_km=[]; stab_gm=[]
for k in SK:
    a_km=[]; a_gm=[]
    for _ in range(25):
        i1=rng.choice(len(Z),int(.8*len(Z)),replace=False); i2=rng.choice(len(Z),int(.8*len(Z)),replace=False)
        sh=np.intersect1d(i1,i2)
        _,C1,_=kmeans(Z[i1],k,seed=1); _,C2,_=kmeans(Z[i2],k,seed=2)
        a_km.append(adjusted_rand_score(((Z[sh][:,None]-C1[None])**2).sum(2).argmin(1),
                                        ((Z[sh][:,None]-C2[None])**2).sum(2).argmin(1)))
        a_gm.append(adjusted_rand_score(GMM(k,seed=1).fit(Z[i1]).predict(Z[sh]),GMM(k,seed=2).fit(Z[i2]).predict(Z[sh])))
    stab_km.append(np.mean(a_km)); stab_gm.append(np.mean(a_gm))
kstab=SK[int(np.argmax(stab_gm))]

fig,ax=plt.subplots(1,3,figsize=(15,4))
ax[0].plot(KS,inertias,"o-",color=BLUE,lw=2); ax[0].axvline(3,color=RED,ls="--",label="elbow ~ 3")
ax[0].set_xlabel("k"); ax[0].set_ylabel("k-means inertia (SSE)"); ax[0].set_title("Elbow method"); ax[0].legend()
ax[1].plot(KS,bics,"o-",color=GREEN,lw=2); ax[1].axvline(kbic,color=RED,ls="--",label=f"BIC minimum: k={kbic}")
ax[1].set_xlabel("k"); ax[1].set_ylabel("GMM BIC (lower=better)"); ax[1].set_title("BIC has an interior minimum"); ax[1].legend()
ax[2].plot(SK,stab_km,"o-",color=BLUE,lw=2,label="k-means"); ax[2].plot(SK,stab_gm,"o-",color=ORANGE,lw=2,label="GMM")
ax[2].axhline(0,color=GREY,lw=.8); ax[2].set_xlabel("k"); ax[2].set_ylabel("agreement across subsamples (ARI)")
ax[2].set_title("Stability: does the partition reproduce?"); ax[2].legend(fontsize=8); ax[2].set_ylim(0,1)
plt.tight_layout(); plt.show()

print(f"The elbow is at 3: inertia falls {1-inertias[2]/inertias[1]:.0%} going to k=3 and only {1-inertias[3]/inertias[2]:.0%} going to k=4.")
print(f"BIC disagrees -- it keeps falling to k={kbic} (BIC {min(bics):.0f}) before turning up. Two defensible criteria, two")
print("different answers, and nothing so far to break the tie: both are in-sample fit measures.")
print("\nStability breaks it. A partition that means something should reappear when the data is resampled:")
for k,a,b in zip(SK,stab_km,stab_gm):
    mark=" <-- elbow" if k==3 else (f" <-- BIC" if k==kbic else "")
    print(f"   k={k}:  k-means ARI {a:.3f}   GMM ARI {b:.3f}{mark}")
print(f"Stability does not name a single k by itself -- k=3, 4 and 5 all reproduce well (ARI {min(min(stab_km[1:4]),min(stab_gm[1:4])):.2f} and above).")
print(f"k=2 is the odd case: stable for the GMM ({stab_gm[0]:.2f}) but not for k-means ({stab_km[0]:.2f}), because with only two")
print("clusters to work with, the direction of the split flips between one resample and the next.")
print(f"What it does decisively is rule out the top of the range: BIC's k={kbic} scores {stab_km[SK.index(kbic)]:.2f} / {stab_gm[SK.index(kbic)]:.2f}, and from k=6 the GMM")
print(f"falls to {min(stab_gm[4:]):.2f} -- barely better than chance agreement. Those extra components land somewhere different")
print("every time they are fitted, which is what an over-split looks like.")
print("So the criteria compose rather than compete: stability eliminates everything above five, and within what")
print("survives the elbow picks three. Neither test alone would have been enough, and BIC alone would have been wrong.")

print(f"\nWhy does BIC want more components than the structure supports? Because a Gaussian mixture will spend components")
print("on non-Gaussian SHAPE, not just on genuinely distinct groups. A control makes the point: draw i.i.d. data from a")
print("single fat-tailed distribution -- one population, no regimes whatsoever -- and score it the same way.")
_r=np.random.default_rng(7); _S=np.cov(Z.T)
for _nu,_tag in [(4.0,"one bivariate t, df=4 (fat tails)"),(np.inf,"one bivariate normal")]:
    _Y=_r.multivariate_normal(np.zeros(2),_S,len(Z))
    if np.isfinite(_nu): _Y=_Y/np.sqrt(_r.chisquare(_nu,len(Z))/_nu)[:,None]
    _Y=StandardScaler().fit_transform(_Y); _b=[GMM(k,seed=0).fit(_Y).bic(_Y) for k in range(1,8)]
    print(f"   {_tag:34s} BIC-minimizing k = {int(np.argmin(_b))+1}")
print("The normal control correctly wants one component. The fat-tailed one -- still a single population, with no")
print("regimes in it at all -- asks for as many components as we ended up using on the real market data, purely to")
print("fit its tails. Daily equity returns are famously fat-tailed, so some of BIC's appetite here is kurtosis being")
print("fitted rather than regimes being found. That is the reason the choice of k rests on the resampling check.")
print("\nThe Bayesian alternative is to stop choosing altogether: the Dirichlet-process mixture (BNP arc) puts a prior")
print("over the number of clusters and INFERS it with uncertainty, rather than fitting a grid and scoring it.")

The elbow is at 3: inertia falls 36% going to k=3 and only 23% going to k=4.
BIC disagrees -- it keeps falling to k=7 (BIC 17224) before turning up. Two defensible criteria, two
different answers, and nothing so far to break the tie: both are in-sample fit measures.

Stability breaks it. A partition that means something should reappear when the data is resampled:
   k=2:  k-means ARI 0.650   GMM ARI 0.950
   k=3:  k-means ARI 0.930   GMM ARI 0.889 <-- elbow
   k=4:  k-means ARI 0.941   GMM ARI 0.848
   k=5:  k-means ARI 0.846   GMM ARI 0.806
   k=6:  k-means ARI 0.684   GMM ARI 0.446
   k=7:  k-means ARI 0.786   GMM ARI 0.571 <-- BIC
   k=8:  k-means ARI 0.667   GMM ARI 0.465
Stability does not name a single k by itself -- k=3, 4 and 5 all reproduce well (ARI 0.81 and above).
k=2 is the odd case: stable for the GMM (0.95) but not for k-means (0.65), because with only two
clusters to work with, the direction of the split flips between one resample and the next.
What it does decisively is rule out the top of the range: BIC's k=7 scores 0.79 / 0.57, and from k=6 the GMM
falls to 0.45 -- barely better than chance agreement. Those extra components land somewhere different
every time they are fitted, which is what an over-split looks like.
So the criteria compose rather than compete: stability eliminates everything above five, and within what
survives the elbow picks three. Neither test alone would have been enough, and BIC alone would have been wrong.

Why does BIC want more components than the structure supports? Because a Gaussian mixture will spend components
on non-Gaussian SHAPE, not just on genuinely distinct groups. A control makes the point: draw i.i.d. data from a
single fat-tailed distribution -- one population, no regimes whatsoever -- and score it the same way.

   one bivariate t, df=4 (fat tails)  BIC-minimizing k = 3

   one bivariate normal               BIC-minimizing k = 1
The normal control correctly wants one component. The fat-tailed one -- still a single population, with no
regimes in it at all -- asks for as many components as we ended up using on the real market data, purely to
fit its tails. Daily equity returns are famously fat-tailed, so some of BIC's appetite here is kurtosis being
fitted rather than regimes being found. That is the reason the choice of k rests on the resampling check.

The Bayesian alternative is to stop choosing altogether: the Dirichlet-process mixture (BNP arc) puts a prior
over the number of clusters and INFERS it with uncertainty, rather than fitting a grid and scoring it.

calm=int(np.argmin(g.mu[:,1]))                                 # calm = the low-volatility cluster
rest=[j for j in range(3) if j!=calm]
sell,rebd=sorted(rest,key=lambda j:ret[lab_g==j].mean())       # of the two turbulent ones, by mean return
order=[calm,sell,rebd]; remap={old:new for new,old in enumerate(order)}
reg=np.array([remap[l] for l in lab_g]); labels=["calm","turbulent sell-off","turbulent rebound"]
print("GMM market regimes (original units):")
for j in range(3):
    m=reg==j
    print(f"  {labels[j]:19s}: mean return {ret[m].mean():+.2f}%,  mean log-vol {lvol[m].mean():+.2f},  "
          f"{m.mean():.0%} of days,  {np.mean(ret[m]<0):.0%} of them down days")
print(f"\nThe two turbulent clusters sit at essentially the same volatility ({lvol[reg==1].mean():+.2f} vs {lvol[reg==2].mean():+.2f} log-vol);")
print(f"what separates them is direction -- {np.mean(ret[reg==1]<0):.0%} of the first are down days against {np.mean(ret[reg==2]<0):.0%} of the second.")
terciles=np.digitize(lvol,np.quantile(lvol,[1/3,2/3]))
print(f"And the partition is not volatility terciles wearing a disguise: its agreement with a plain three-way split of")
print(f"log-volatility is only ARI {adjusted_rand_score(reg,terciles):.3f}. The return axis is doing real work.")

fig,ax=plt.subplots(figsize=(12,4)); cols=[GREEN,ORANGE,RED]
av=np.sqrt(d["rv"].values)*100*np.sqrt(252)
for j in range(3): m=reg==j; ax.scatter(dates[m],av[m],s=5,color=cols[j],label=labels[j])
ax.set_ylabel("annualized realized vol (%)"); ax.set_title("Market regimes over time (colored by GMM cluster)"); ax.legend(markerscale=2,fontsize=8)
plt.tight_layout(); plt.show()

turb=reg>0; yr=dates.dt.year.values
print(f"Turbulent share by year (unconditional {turb.mean():.0%}):")
print("   "+"  ".join(f"{y}:{turb[yr==y].mean():>4.0%}" for y in range(2000,2007)))
print("   "+"  ".join(f"{y}:{turb[yr==y].mean():>4.0%}" for y in range(2007,2014)))
_top=sorted(range(2000,2014),key=lambda y:-turb[yr==y].mean())[:3]
print(f"The turbulent regimes concentrate in two distinct episodes, not one: {_top[0]} ({turb[yr==_top[0]].mean():.0%} of days),")
print(f"{_top[1]} ({turb[yr==_top[1]].mean():.0%}) and {_top[2]} ({turb[yr==_top[2]].mean():.0%}) -- the financial crisis and its aftermath, and the dot-com")
print(f"unwind, which is fully comparable rather than a footnote to it. The calm years {sorted(range(2000,2014),key=lambda y:turb[yr==y].mean())[0]}-2006 sit")
print(f"near {min(turb[yr==y].mean() for y in [2004,2005,2006]):.0%}. Discovered with no dates or labels supplied -- the unsupervised counterpart to the")
print("crisis dating a human would do by hand.")

GMM market regimes (original units):
  calm               : mean return +0.08%,  mean log-vol -0.46,  66% of days,  43% of them down days
  turbulent sell-off : mean return -1.06%,  mean log-vol +0.31,  24% of days,  76% of them down days
  turbulent rebound  : mean return +2.06%,  mean log-vol +0.47,  10% of days,  1% of them down days

The two turbulent clusters sit at essentially the same volatility (+0.31 vs +0.47 log-vol);
what separates them is direction -- 76% of the first are down days against 1% of the second.
And the partition is not volatility terciles wearing a disguise: its agreement with a plain three-way split of
log-volatility is only ARI 0.270. The return axis is doing real work.

Turbulent share by year (unconditional 34%):
   2000: 56%  2001: 46%  2002: 64%  2003: 37%  2004:  8%  2005:  4%  2006:  6%
   2007: 27%  2008: 68%  2009: 62%  2010: 32%  2011: 42%  2012: 16%  2013:  9%
The turbulent regimes concentrate in two distinct episodes, not one: 2008 (68% of days),
2002 (64%) and 2009 (62%) -- the financial crisis and its aftermath, and the dot-com
unwind, which is fully comparable rather than a footnote to it. The calm years 2005-2006 sit
near 4%. Discovered with no dates or labels supplied -- the unsupervised counterpart to the
crisis dating a human would do by hand.

Unsupervised Learning I — Clustering: k-means & Gaussian mixtures¶

Built from scratch, validated against scikit-learn, applied to market regimes¶

1. The data — market states without labels¶

2. k-means from scratch¶

3. Gaussian mixtures from scratch (EM)¶

4. How many clusters? — elbow, BIC, and the Bayesian answer¶

5. The regimes, interpreted — and over time¶

6. Summary¶