Kapitel 18: Finanzdaten-Modellierung — Renditen, Kovarianz und Shrinkage

#!/usr/bin/env python3

# Renditen_Vergleich.py
"""
Kapitel Finanzdaten: Diskrete vs. logarithmische Renditen.
Zeigt an einem simulierten Kursverlauf, welche Eigenschaft wo gilt -
und warum das arithmetische Mittel diskreter Renditen in die Irre fuehrt.
"""

import numpy as np
import pandas as pd

if __name__ == "__main__":
    kurse = np.array([100.0, 150.0, 100.0, 120.0, 90.0, 135.0])
    tage = [f"t{i}" for i in range(len(kurse))]

    diskret = kurse[1:] / kurse[:-1] - 1.0
    logarithmisch = np.log(kurse[1:] / kurse[:-1])

    print("=" * 74)
    print("  DISKRETE UND LOGARITHMISCHE RENDITEN IM VERGLEICH")
    print("=" * 74)
    tabelle = pd.DataFrame({
        "Periode": [f"{tage[i]} -> {tage[i+1]}" for i in range(len(diskret))],
        "Kurs von": kurse[:-1],
        "Kurs bis": kurse[1:],
        "diskret R": [f"{r*100:+7.2f} %" for r in diskret],
        "log r": [f"{r:+8.4f}" for r in logarithmisch],
    })
    print(tabelle.to_string(index=False))

    # --- Zeitliche Aggregation ------------------------------------------
    gesamt_wahr = kurse[-1] / kurse[0] - 1.0
    summe_log = logarithmisch.sum()
    aus_log_zurueck = np.exp(summe_log) - 1.0
    summe_diskret = diskret.sum()

    print("\n--- Aggregation ueber die Zeit ---")
    print(f"  Tatsaechliche Gesamtrendite:        {gesamt_wahr*100:+8.2f} %")
    print(f"  Summe der Log-Renditen -> exp()-1:  {aus_log_zurueck*100:+8.2f} %  "
          f"{'KORREKT' if abs(aus_log_zurueck - gesamt_wahr) < 1e-9 else 'falsch'}")
    print(f"  Summe der diskreten Renditen:       {summe_diskret*100:+8.2f} %  FALSCH")

    # --- Mittelwerte ------------------------------------------------------
    arithmetisch = diskret.mean()
    geometrisch = np.prod(1 + diskret) ** (1 / len(diskret)) - 1
    aus_log = np.exp(logarithmisch.mean()) - 1

    print("\n--- Welcher Mittelwert ist der richtige? ---")
    print(f"  Arithmetisches Mittel (diskret):    {arithmetisch*100:+8.2f} %  "
          f"-> ueberschaetzt")
    print(f"  Geometrisches Mittel:               {geometrisch*100:+8.2f} %  -> korrekt")
    print(f"  exp(Mittel der Log-Renditen) - 1:   {aus_log*100:+8.2f} %  "
          f"-> identisch zum geometrischen")

    # Probe: Endkapital mit dem jeweiligen Mittelwert hochgerechnet
    n = len(diskret)
    print(f"\n  Probe - Startkapital 100 EUR ueber {n} Perioden:")
    print(f"    tatsaechlich:                     {kurse[-1]:8.2f} EUR")
    print(f"    mit arithm. Mittel hochgerechnet: {100*(1+arithmetisch)**n:8.2f} EUR")
    print(f"    mit geom.  Mittel hochgerechnet:  {100*(1+geometrisch)**n:8.2f} EUR")

    print("\n--- Annualisierung (252 Handelstage) ---")
    print(f"  Volatilitaet aus Log-Renditen: "
          f"{logarithmisch.std(ddof=1)*np.sqrt(252)*100:.2f} % p.a.")
    print("  (Die Wurzel-Zeit-Regel gilt nur bei unabhaengigen Renditen -")
    print("   fuer reale Maerkte ist sie eine Naeherung.)")
    print("=" * 74)

#!/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)

#!/usr/bin/env python3

# Finanzdaten_Ledoit_Wolf.py
"""
Kapitel Finanzdaten: Automatisierte Finanzdaten-Pipeline und Kovarianz-Shrinkage.

Eigenschaften:
  * Spaltenreihenfolge wird erzwungen und per assert geprueft
  * beide Shrinkage-Ziele werden berechnet und verglichen
  * Datenqualitaetspruefungen (Luecken, Ausreisser, Mindestlaenge)
  * Eigenwertspektrum wird ausgewiesen, nicht nur die Konditionszahl
"""

import numpy as np
import pandas as pd
import yfinance as yf
from sklearn.covariance import LedoitWolf

HANDELSTAGE = 252


def lade_kurse(tickers: list[str], start, ende) -> pd.DataFrame:
    """Laedt bereinigte Schlusskurse und garantiert die Spaltenreihenfolge."""
    print(f"Lade {len(tickers)} Ticker von {start} bis {ende} ...")
    roh = yf.download(tickers, start=start, end=ende, auto_adjust=True, progress=False)
    if roh.empty:
        raise SystemExit("Download fehlgeschlagen (Netz, Ticker oder Rate-Limit pruefen).")

    if isinstance(roh.columns, pd.MultiIndex):
        # [tickers] erzwingt die gewuenschte Spaltenreihenfolge
        kurse = roh["Close"][tickers].dropna()
    else:
        kurse = roh[["Close"]].dropna()
        kurse.columns = tickers

    assert list(kurse.columns) == tickers, (
        f"Spaltenreihenfolge weicht ab!\n  erwartet: {tickers}\n"
        f"  erhalten: {list(kurse.columns)}")
    return kurse


def pruefe_datenqualitaet(kurse: pd.DataFrame, mindest_tage: int = 200) -> None:
    """Faengt die haeufigsten Datenprobleme ab, bevor sie ins Modell gelangen."""
    T, N = kurse.shape
    print(f"\nDatenqualitaet: {T} Handelstage, {N} Titel (T/N = {T/N:.1f})")
    if T < mindest_tage:
        raise SystemExit(f"Zu wenige Beobachtungen ({T} < {mindest_tage}).")
    if T < N:
        print("  WARNUNG: T < N - die Stichprobenkovarianz ist singulaer!")
    elif T < 3 * N:
        print("  WARNUNG: T < 3N - erhebliches Schaetzrauschen zu erwarten.")

    renditen = kurse.pct_change().dropna()
    extreme = (renditen.abs() > 0.25).sum().sum()
    if extreme:
        print(f"  Hinweis: {extreme} Tagesrenditen ueber 25 % "
              f"(Splits, Sondersituationen oder Datenfehler pruefen).")
    luecken = kurse.isna().sum().sum()
    print(f"  Fehlende Werte nach dropna: {luecken}")


def ziel_konstante_korrelation(renditen: np.ndarray) -> np.ndarray:
    """
    Shrinkage-Ziel nach Ledoit/Wolf 2003: Einzelvarianzen behalten,
    alle Korrelationen durch ihren Mittelwert ersetzen.
    (sklearn nutzt stattdessen die skalierte Einheitsmatrix - siehe oben.)
    """
    S = np.cov(renditen, rowvar=False, ddof=1)
    d = np.sqrt(np.diag(S))
    korrelation = S / np.outer(d, d)
    n = len(S)
    r_quer = (korrelation.sum() - n) / (n * (n - 1))       # Mittel ohne Diagonale
    F = r_quer * np.outer(d, d)
    np.fill_diagonal(F, np.diag(S))
    return F


def analysiere(kurse: pd.DataFrame):
    renditen = kurse.pct_change().dropna()
    T, N = renditen.shape

    mu = renditen.mean().values * HANDELSTAGE
    sigma_stichprobe = renditen.cov().values * HANDELSTAGE

    lw = LedoitWolf(assume_centered=False).fit(renditen.values)
    sigma_lw = lw.covariance_ * HANDELSTAGE
    delta = lw.shrinkage_

    print("\n" + "=" * 84)
    print("         FINANZDATEN-PIPELINE UND MATRIX-KONDITIONIERUNG")
    print("=" * 84)
    print(f"Beobachtungen T:            {T}")
    print(f"Titel N:                    {N}")
    print(f"Optimales Shrinkage delta:  {delta:.4f} "
          f"({delta*100:.1f} % Gewicht auf dem strukturierten Ziel)")

    eig_stich = np.linalg.eigvalsh(sigma_stichprobe)
    eig_lw = np.linalg.eigvalsh(sigma_lw)

    print("\n--- Eigenwertspektrum (annualisiert) ---")
    print(f"{'':<14} {'kleinster':>12} {'Median':>12} {'groesster':>12} "
          f"{'Kondition':>12}")
    for name, eig in [("Stichprobe", eig_stich), ("Ledoit-Wolf", eig_lw)]:
        kondition = eig.max() / max(eig.min(), 1e-12)
        print(f"{name:<14} {eig.min():>12.6f} {np.median(eig):>12.6f} "
              f"{eig.max():>12.6f} {kondition:>12.1f}")

    print("\nDie Konditionszahl misst, wie stark sich kleine Datenaenderungen auf")
    print("die Inverse auswirken. Je kleiner, desto stabiler die Portfoliogewichte.")

    # --- Vergleich der beiden Shrinkage-Ziele ----------------------------
    S_taeglich = np.cov(renditen.values, rowvar=False, ddof=1)
    F_identitaet = np.eye(N) * np.trace(S_taeglich) / N
    F_korrelation = ziel_konstante_korrelation(renditen.values)

    print("\n--- Die beiden Shrinkage-Ziele im Vergleich ---")
    for name, F in [("skalierte Einheitsmatrix (sklearn)", F_identitaet),
                    ("konstante Korrelation (LW 2003)", F_korrelation)]:
        gemischt = ((1 - delta) * S_taeglich + delta * F) * HANDELSTAGE
        eig = np.linalg.eigvalsh(gemischt)
        print(f"  {name:<36}: Kondition {eig.max()/max(eig.min(),1e-12):8.1f}, "
              f"kleinster EW {eig.min():.6f}")

    print("\n--- Erwartete Renditen (annualisiert) ---")
    uebersicht = pd.DataFrame({
        "Ticker": kurse.columns,
        "Rendite p.a.": [f"{r*100:+7.2f} %" for r in mu],
        "Volatilitaet p.a.": [f"{np.sqrt(sigma_lw[i, i])*100:6.2f} %" for i in range(N)],
    })
    print(uebersicht.to_string(index=False))
    print("=" * 84)
    return renditen, mu, sigma_stichprobe, sigma_lw


if __name__ == "__main__":
    UNIVERSUM = ["AAPL", "MSFT", "NVDA", "AMZN",   # Technologie
                 "JNJ", "PFE",                      # Gesundheit
                 "JPM", "GS",                       # Banken
                 "XOM", "CVX"]                      # Energie

    # Relatives 2-Jahres-Fenster, damit das Beispiel nicht "altert"
    ende = pd.Timestamp.today().normalize()
    start = ende - pd.DateOffset(years=2)

    kurse = lade_kurse(UNIVERSUM, start.strftime("%Y-%m-%d"), ende.strftime("%Y-%m-%d"))
    pruefe_datenqualitaet(kurse)
    analysiere(kurse)

#!/usr/bin/env python3

# Kovarianz_Falle.py
"""
Kapitel Finanzdaten: Warum die historische Kovarianzmatrix den Optimierer belaeugt.

Ein Portfolio aus 30 Anlagen soll das Risiko minimieren. Geschaetzt wird die
Kovarianzmatrix aus 20 Beobachtungen - in der Praxis ein voellig ueblicher
Fall (30 Titel, ein Quartal Wochendaten).

Das Ergebnis ist ein Portfolio mit einem gemessenen Risiko von 0,000 %. Der
Optimierer haelt es fuer risikofrei. Ausserhalb des Schaetzzeitraums ist es
das schlechteste der drei untersuchten Portfolios - schlechter sogar als
blosse Gleichgewichtung.

Der Grund ist nicht Zufall, sondern Struktur: Bei N Anlagen und T < N
Beobachtungen hat die Stichproben-Kovarianzmatrix hoechstens Rang T-1. Es gibt
dann ganze Richtungen im Gewichtsraum, in denen sie exakt null Varianz misst -
Richtungen, die es in Wirklichkeit nicht gibt. Der Optimierer findet sie
zuverlaessig, denn er sucht ja genau danach.

Die Daten sind synthetisch (ein Marktfaktor plus Eigenrauschen), damit das
Ergebnis reproduzierbar ist und kein Internetzugang noetig wird. Mit echten
Kursdaten sieht es genauso aus.

Benoetigt: numpy, cvxpy, scikit-learn
"""

from __future__ import annotations

import numpy as np
import cvxpy as cp
from sklearn.covariance import LedoitWolf

ANZAHL_ANLAGEN = 30
BEOBACHTUNGEN_SCHAETZUNG = 20        # weniger als Anlagen - genau darum geht es
BEOBACHTUNGEN_TEST = 2000            # der "spaetere Verlauf"

RNG = np.random.default_rng(11)
# Jede Anlage reagiert unterschiedlich stark auf den Gesamtmarkt.
BETA = RNG.uniform(0.6, 1.4, ANZAHL_ANLAGEN)


def erzeuge_renditen(perioden: int) -> np.ndarray:
    """Ein Marktfaktor, auf den alle reagieren, plus Eigenrauschen je Anlage."""
    markt = RNG.normal(0.0, 0.010, perioden)
    eigen = RNG.normal(0.0, 0.012, (perioden, ANZAHL_ANLAGEN))
    return np.outer(markt, BETA) + eigen


def minimales_risiko(kovarianz: np.ndarray) -> tuple[np.ndarray, float]:
    """Minimum-Varianz-Portfolio: Gewichte summieren sich zu 1, sonst frei.

    cp.psd_wrap() sagt CVXPY: 'Ich weiss, dass diese Matrix numerisch nicht
    exakt positiv semidefinit ist - rechne trotzdem.' Das ist hier bewusst so
    gewaehlt, denn genau diese Eigenschaft wollen wir vorfuehren. Im
    Produktivbetrieb waere die DCP-Fehlermeldung das Warnsignal, dem man
    nachgehen muss (siehe den Denkfehler im Kapitel Quadratische und
    nichtlineare Optimierung).
    """
    gewichte = cp.Variable(ANZAHL_ANLAGEN)
    problem = cp.Problem(cp.Minimize(cp.quad_form(gewichte, cp.psd_wrap(kovarianz))),
                         [cp.sum(gewichte) == 1])
    problem.solve()
    return np.array(gewichte.value).ravel(), max(problem.value, 0.0)


if __name__ == "__main__":
    schaetzzeitraum = erzeuge_renditen(BEOBACHTUNGEN_SCHAETZUNG)
    spaeterer_verlauf = erzeuge_renditen(BEOBACHTUNGEN_TEST)

    stichprobe = np.cov(schaetzzeitraum, rowvar=False)
    ledoit_wolf = LedoitWolf().fit(schaetzzeitraum)

    print("=" * 78)
    print("  DIE KOVARIANZ-FALLE")
    print("=" * 78)
    print(f"{ANZAHL_ANLAGEN} Anlagen, geschaetzt aus {BEOBACHTUNGEN_SCHAETZUNG} "
          f"Beobachtungen.")
    print(f"Rang der Stichprobenmatrix: "
          f"{np.linalg.matrix_rank(stichprobe)} - noetig waeren {ANZAHL_ANLAGEN}.")
    print(f"Kleinster Eigenwert: {np.linalg.eigvalsh(stichprobe).min():.2e}")
    print()
    print("Es gibt also Richtungen, in denen diese Matrix EXAKT null Varianz")
    print("misst. Genau dort sucht ein Minimum-Varianz-Optimierer.")
    print()

    print(f"{'Schaetzer':<26} {'Risiko im Zeitraum':>19} "
          f"{'Risiko spaeter':>16} {'Hebel':>7}")
    print("-" * 78)

    ergebnisse = {}
    for name, matrix in [
            ("Stichprobe (roh)", stichprobe),
            (f"Ledoit-Wolf (d={ledoit_wolf.shrinkage_:.2f})", ledoit_wolf.covariance_)]:
        gewichte, varianz_intern = minimales_risiko(matrix)
        risiko_spaeter = float(np.std(spaeterer_verlauf @ gewichte))
        hebel = float(np.abs(gewichte).sum())
        ergebnisse[name] = (risiko_spaeter, gewichte)
        print(f"{name:<26} {np.sqrt(varianz_intern) * 100:>18.3f} % "
              f"{risiko_spaeter * 100:>15.3f} % {hebel:>7.1f}")

    gleich = np.full(ANZAHL_ANLAGEN, 1.0 / ANZAHL_ANLAGEN)
    print(f"{'Gleichgewichtung 1/N':<26} "
          f"{float(np.std(schaetzzeitraum @ gleich)) * 100:>18.3f} % "
          f"{float(np.std(spaeterer_verlauf @ gleich)) * 100:>15.3f} % {1.0:>7.1f}")
    print("-" * 78)

    print("\nDrei Beobachtungen:")
    print("1. Die rohe Schaetzung meldet 0,000 % Risiko - und liefert spaeter das")
    print("   SCHLECHTESTE Ergebnis der drei. Sie hat kein Portfolio optimiert,")
    print("   sondern eine Luecke in den eigenen Daten gefunden.")
    print("2. Die Gleichgewichtung, die gar nicht optimiert, schlaegt den rohen")
    print("   Optimierer deutlich. Das ist kein Zufallsbefund, sondern in der")
    print("   Literatur breit belegt.")
    print("3. Ledoit-Wolf ist im Schaetzzeitraum EHRLICHER (0,441 statt 0,000 %)")
    print("   und spaeter am besten. Der Hebel faellt von 3,8 auf 1,8 - das")
    print("   Portfolio wird auch praktisch handelbarer.")
    print()
    print("Merksatz: Ein Optimierer glaubt seinen Eingabedaten vollstaendig.")
    print("Je mehr Freiheit Sie ihm geben, desto gruendlicher findet er deren")
    print("Fehler. Das gilt fuer Kovarianzmatrizen wie fuer Lieferzeiten,")
    print("Ausfallraten und jede andere geschaetzte Groesse.")
    print("=" * 78)