#!/usr/bin/env python3

# Schaetzrauschen_Demo.py
"""
Kapitel Finanzdaten: Der Error-Maximizer-Effekt im Experiment.

Aufbau: Wir KENNEN die wahre Kovarianzmatrix, weil wir die Daten selbst
erzeugen. Dann schaetzen wir sie aus endlich vielen Beobachtungen und
vergleichen, wie gut das daraus optimierte Portfolio in der WAHREN Welt
abschneidet - mit und ohne Shrinkage.

Das ist der einzige saubere Weg, Schaetzfehler zu messen: Man braucht eine
Wahrheit zum Vergleich, und die gibt es nur in der Simulation.
"""

import numpy as np
from sklearn.covariance import LedoitWolf


def erzeuge_wahre_kovarianz(n, rng):
    """Ein-Faktor-Modell: gemeinsamer Marktfaktor plus titelspezifisches Rauschen."""
    beta = rng.uniform(0.6, 1.4, size=n)             # Marktsensitivitaeten
    var_markt = 0.04                                  # Marktvarianz p.a.
    var_spezifisch = rng.uniform(0.01, 0.09, size=n)  # idiosynkratische Varianz
    return np.outer(beta, beta) * var_markt + np.diag(var_spezifisch)


def minimum_varianz_gewichte(sigma):
    """Analytische GMV-Loesung: w = Sigma^-1 1 / (1' Sigma^-1 1). Ohne Restriktionen."""
    eins = np.ones(len(sigma))
    try:
        loesung = np.linalg.solve(sigma, eins)
    except np.linalg.LinAlgError:
        loesung = np.linalg.pinv(sigma) @ eins        # falls singulaer
    return loesung / loesung.sum()


def portfoliovolatilitaet(w, sigma):
    return float(np.sqrt(w @ sigma @ w))


if __name__ == "__main__":
    rng = np.random.default_rng(2026)
    N = 40                                            # Anzahl Titel
    WIEDERHOLUNGEN = 200

    sigma_wahr = erzeuge_wahre_kovarianz(N, rng)
    w_ideal = minimum_varianz_gewichte(sigma_wahr)
    vola_ideal = portfoliovolatilitaet(w_ideal, sigma_wahr)

    print("=" * 92)
    print("  ERROR-MAXIMIZER: WIE TEUER IST SCHAETZRAUSCHEN?")
    print("=" * 92)
    print(f"Aufbau: {N} Titel, wahre Kovarianz bekannt (Ein-Faktor-Modell).")
    print(f"Bestmoegliche Volatilitaet bei perfektem Wissen: "
          f"{vola_ideal*100:.2f} % p.a.\n")

    print(f"{'T (Tage)':>9} | {'T/N':>5} | {'Stichprobe':>22} | {'Ledoit-Wolf':>22} | "
          f"{'Shrink':>7}")
    print(f"{'':>9} | {'':>5} | {'Vola     Aufschlag':>22} | "
          f"{'Vola     Aufschlag':>22} | {'delta':>7}")
    print("-" * 92)

    for T in [60, 120, 252, 504, 1260]:
        vola_stichprobe, vola_lw, deltas = [], [], []

        for _ in range(WIEDERHOLUNGEN):
            # Daten aus der WAHREN Verteilung ziehen (taeglich)
            daten = rng.multivariate_normal(np.zeros(N), sigma_wahr / 252, size=T)

            # (a) Stichproben-Kovarianz
            s_stich = np.cov(daten, rowvar=False, ddof=1) * 252
            w_stich = minimum_varianz_gewichte(s_stich)

            # (b) Ledoit-Wolf-Shrinkage
            lw = LedoitWolf().fit(daten)
            s_lw = lw.covariance_ * 252
            w_lw = minimum_varianz_gewichte(s_lw)
            deltas.append(lw.shrinkage_)

            # Bewertung IMMER mit der wahren Kovarianz - das ist der Punkt
            vola_stichprobe.append(portfoliovolatilitaet(w_stich, sigma_wahr))
            vola_lw.append(portfoliovolatilitaet(w_lw, sigma_wahr))

        m_stich, m_lw = np.mean(vola_stichprobe), np.mean(vola_lw)
        print(f"{T:>9} | {T/N:>5.1f} | {m_stich*100:>8.2f} % "
              f"{(m_stich/vola_ideal-1)*100:>+9.1f} % | "
              f"{m_lw*100:>8.2f} % {(m_lw/vola_ideal-1)*100:>+9.1f} % | "
              f"{np.mean(deltas):>7.3f}")

    print("-" * 92)
    print("'Aufschlag' = wie viel mehr Risiko das Portfolio in der WAHREN Welt traegt,")
    print("verglichen mit dem Portfolio bei perfektem Wissen.")
    print("Lesart: Je kleiner T/N, desto teurer das Schaetzrauschen - und desto mehr")
    print("schrumpft Ledoit-Wolf (delta steigt).")
    print("=" * 92)