L08 — Componenti principali e numero effettivo di scommesse

import datetime as dt

import matplotlib.pyplot as plt
import numpy as np
import polars as pl

from cvbook.dati import carica
from cvbook.metriche import rendimenti


def allinea(nomi: list[str], da: dt.date = dt.date(2020, 9, 1)) -> np.ndarray:
    """Rendimenti giornalieri delle serie, sui soli giorni presenti in tutte."""
    pezzi = [
        carica(n).filter(pl.col("data") >= da).select(["data", "chiusura"]).rename({"chiusura": n})
        for n in nomi
    ]
    base = pezzi[0]
    for p in pezzi[1:]:
        base = base.join(p, on="data")
    base = base.sort("data")
    return np.column_stack([rendimenti(base[n].to_numpy()) for n in nomi])

NOMI = ["btcusdt", "ethusdt", "solusdt"]
# PROVA / TRY: togli "solusdt" (esercizio 1) · aggiungi "ftsemib"/"eni"
# (aggiungili anche a avvio.prepara([...]))

M = allinea(NOMI)
C = np.corrcoef(M.T)
autovalori = np.linalg.eigvalsh(C)[::-1]
quote = autovalori / autovalori.sum()
cumulata = np.cumsum(quote)

with avvio.figura("schermo"):
    fig, (sx, dx) = plt.subplots(1, 2, figsize=(11, 4))
    x = np.arange(1, len(quote) + 1)
    sx.bar(x, quote * 100)
    for k, q in enumerate(quote):
        sx.annotate(f"{q:.0%}", xy=(k + 1, q * 100), xytext=(0, 4),
                    textcoords="offset points", ha="center")
    sx.set_xticks(x)
    sx.set_xlabel("Componente")
    sx.set_ylabel("Varianza spiegata (%)")

    dx.plot(x, cumulata * 100, marker="o")
    dx.axhline(90, linestyle=":", linewidth=1.2)
    dx.set_xticks(x)
    dx.set_xlabel("Componenti usate")
    dx.set_ylabel("Varianza cumulata (%)")
    plt.show()

for k, (q, c) in enumerate(zip(quote, cumulata), start=1):
    print(f"componente {k}: {q:6.1%}   cumulata {c:6.1%}")
print(f"\ncomponenti necessarie per arrivare al 90%: {int(np.searchsorted(cumulata, 0.90)) + 1}")
print(f"giorni usati: {len(M)} — il calendario delle sole cripto, che non chiudono mai")
print("Il capitolo stampa 78,7% perche' misura queste stesse tre serie sui soli "
      "giorni in cui e' aperta anche la borsa di Milano: e' l'unico modo di "
      "confrontarle con il paniere esteso. Due decimi di punto di differenza, "
      "e il numero effettivo di scommesse — 1,55 — e' identico.")

def numero_effettivo(nomi: list[str]) -> tuple[float, float]:
    m = allinea(nomi)
    c = np.corrcoef(m.T)
    v = np.linalg.eigvalsh(c)[::-1]
    q = v / v.sum()
    rho = float(c[np.triu_indices(len(nomi), 1)].mean())
    da_correlazione = 1.0 / (1 / len(nomi) + (1 - 1 / len(nomi)) * rho)
    return float(1.0 / np.sum(q**2)), da_correlazione


da_componenti, da_correlazione = numero_effettivo(NOMI)
print(f"numero effettivo di scommesse (componenti):        {da_componenti:.2f}")
print(f"numero effettivo di scommesse (correlazione media): {da_correlazione:.2f}")
print("\nDue metodi con assunzioni diverse. Il secondo assume che tutte le coppie "
      "abbiano la stessa correlazione; il primo vede la struttura reale. Quando "
      "concordano sull'ordine di grandezza, la conclusione e' molto piu' solida.")

print(f"{'portafoglio':>34s} {'asset':>6s} {'scommesse effettive':>20s}")
for k in range(1, len(NOMI) + 1):
    sottoinsieme = NOMI[:k]
    if k == 1:
        effettivo = 1.0
    else:
        effettivo, _ = numero_effettivo(sottoinsieme)
    print(f"{' + '.join(sottoinsieme):>34s} {k:6d} {effettivo:20.2f}")

print("\nAggiungere asset dello stesso tipo non aggiunge dimensioni: aggiunge "
      "costi e cose da seguire.")

FINESTRA = 250  # PROVA / TRY: 60 (vedi esercizio 2) · 250
prime = []
for i in range(FINESTRA, len(M), 25):
    blocco = M[i - FINESTRA:i]
    v = np.linalg.eigvalsh(np.corrcoef(blocco.T))[::-1]
    prime.append(v[0] / v.sum())

prime = np.array(prime)
with avvio.figura("schermo"):
    fig, ax = plt.subplots(figsize=(9, 3.5))
    ax.plot(prime * 100, marker="o", markersize=3)
    ax.set_xlabel(f"Finestre mobili di {FINESTRA} giorni")
    ax.set_ylabel("Varianza spiegata dalla prima componente (%)")
    ax.set_ylim(0, 100)
    plt.show()

print(f"prima componente: minimo {prime.min():.1%}, massimo {prime.max():.1%}, "
      f"media {prime.mean():.1%}")