#!/usr/bin/env python3

# Renditeschaetzung_Falle.py
"""
Kapitel Markowitz: Warum geschaetzte Renditen noch gefaehrlicher sind als geschaetzte
Kovarianzen.

Das Kapitel Finanzdaten hat gezeigt, was eine schlecht geschaetzte Kovarianzmatrix
anrichtet. Bei den erwarteten RENDITEN ist es schlimmer - aus einem einfachen
statistischen Grund: Um eine Rendite von 0,04 % pro Tag von null zu
unterscheiden, braucht man bei 1,2 % Tagesschwankung ueber tausend
Beobachtungen. Kovarianzen konvergieren dagegen deutlich schneller.

Das Experiment ist bewusst so gebaut, dass die Antwort feststeht: ALLE zwoelf
Anlagen haben exakt dieselbe wahre Rendite und dieselbe Schwankung. Ein
korrektes Ergebnis waere also Gleichgewichtung. Was der Optimierer stattdessen
tut, zeigt der Lauf.

Wiederholt wird 200 mal, damit das Ergebnis nicht von einer Zufallsstichprobe
abhaengt.

Benoetigt: numpy, cvxpy
"""

from __future__ import annotations

import numpy as np
import cvxpy as cp

ANZAHL_ANLAGEN = 12
BEOBACHTUNGEN = 250              # ein Jahr Tagesdaten
BEOBACHTUNGEN_SPAETER = 2000     # der "spaetere Verlauf"
WIEDERHOLUNGEN = 200

# Die Wahrheit, die der Optimierer nicht kennt: alle Anlagen sind gleich.
WAHRE_RENDITE = 0.0004           # 0,04 % pro Tag
WAHRE_SCHWANKUNG = 0.012         # 1,2 % pro Tag
RISIKOAVERSION = 10.0

RNG = np.random.default_rng(5)


def baue_optimierer() -> tuple[cp.Problem, cp.Variable, cp.Parameter, cp.Parameter]:
    """Baut das Markowitz-Problem EINMAL mit Parametern.

    Bei 200 Wiederholungen ist das der Unterschied zwischen Sekunden und
    Minuten: cp.Parameter erlaubt es, nur die Daten zu tauschen, statt den
    Ausdrucksbaum jedes Mal neu zu kompilieren (siehe Kapitel Oekosystem).
    """
    gewichte = cp.Variable(ANZAHL_ANLAGEN, nonneg=True)
    renditen = cp.Parameter(ANZAHL_ANLAGEN)
    kovarianz = cp.Parameter((ANZAHL_ANLAGEN, ANZAHL_ANLAGEN), PSD=True)
    problem = cp.Problem(
        cp.Maximize(renditen @ gewichte
                    - RISIKOAVERSION * cp.quad_form(gewichte, kovarianz)),
        [cp.sum(gewichte) == 1])
    return problem, gewichte, renditen, kovarianz


def erzeuge(perioden: int) -> np.ndarray:
    """Renditen, bei denen alle Anlagen identisch verteilt und unabhaengig sind."""
    return RNG.normal(WAHRE_RENDITE, WAHRE_SCHWANKUNG,
                      (perioden, ANZAHL_ANLAGEN))


if __name__ == "__main__":
    problem, gewichte, renditen, kovarianz = baue_optimierer()

    im_zeitraum, spaeter, spaeter_gleich = [], [], []
    groesstes_gewicht, schaetzspanne = [], []

    for _ in range(WIEDERHOLUNGEN):
        schaetzdaten = erzeuge(BEOBACHTUNGEN)
        spaetere_daten = erzeuge(BEOBACHTUNGEN_SPAETER)

        geschaetzte_rendite = schaetzdaten.mean(axis=0)
        geschaetzte_kovarianz = np.cov(schaetzdaten, rowvar=False)

        renditen.value = geschaetzte_rendite
        # Symmetrisieren und minimal anheben: numerisches Rauschen kann die
        # Matrix sonst knapp unter die PSD-Grenze druecken (Kapitel Fundament).
        kovarianz.value = ((geschaetzte_kovarianz + geschaetzte_kovarianz.T) / 2
                           + 1e-10 * np.eye(ANZAHL_ANLAGEN))
        problem.solve()

        w = np.array(gewichte.value).ravel()
        im_zeitraum.append(float((schaetzdaten @ w).mean()))
        spaeter.append(float((spaetere_daten @ w).mean()))
        spaeter_gleich.append(float(spaetere_daten.mean()))
        groesstes_gewicht.append(float(w.max()))
        schaetzspanne.append(float(geschaetzte_rendite.max()
                                   - geschaetzte_rendite.min()))

    prozent = lambda werte: float(np.mean(werte)) * 100

    print("=" * 78)
    print("  WENN DER OPTIMIERER GESCHAETZTE RENDITEN GLAUBT")
    print("=" * 78)
    print(f"{ANZAHL_ANLAGEN} Anlagen, {WIEDERHOLUNGEN} Wiederholungen, "
          f"je {BEOBACHTUNGEN} Beobachtungen zur Schaetzung.")
    print(f"In Wahrheit haben ALLE dieselbe Tagesrendite von "
          f"{WAHRE_RENDITE * 100:.4f} % und dieselbe Schwankung.")
    print("Die richtige Antwort waere also: gleichgewichten.")
    print()
    print(f"Mittlere Spanne der geschaetzten Renditen: "
          f"{prozent(schaetzspanne):.4f} Prozentpunkte")
    print(f"  -> Das ist das {prozent(schaetzspanne) / (WAHRE_RENDITE * 100):.1f}-Fache "
          f"des wahren Werts. Aus reinem Rauschen.")
    print(f"Mittleres groesstes Einzelgewicht: "
          f"{prozent(groesstes_gewicht):.1f} % "
          f"(bei Gleichgewichtung waeren es {100 / ANZAHL_ANLAGEN:.1f} %)")
    print()
    print(f"{'':<38} {'Tagesrendite':>14}")
    print("-" * 78)
    print(f"{'Wahrheit (alle Anlagen)':<38} {WAHRE_RENDITE * 100:>13.4f} %")
    print(f"{'Optimierer IM Schaetzzeitraum':<38} {prozent(im_zeitraum):>13.4f} %"
          f"   <- Illusion")
    print(f"{'Optimierer AUSSERHALB':<38} {prozent(spaeter):>13.4f} %")
    print(f"{'Gleichgewichtung AUSSERHALB':<38} {prozent(spaeter_gleich):>13.4f} %")
    print("-" * 78)

    faktor = prozent(im_zeitraum) / (WAHRE_RENDITE * 100)
    vorsprung = prozent(spaeter) - prozent(spaeter_gleich)
    print(f"\nIm Schaetzzeitraum verspricht der Optimierer das "
          f"{faktor:.1f}-Fache der wahren Rendite.")
    print(f"Ausserhalb bleibt davon nichts: Er liegt {abs(vorsprung):.4f} "
          f"Prozentpunkte")
    print(f"{'schlechter' if vorsprung < 0 else 'besser'} als blosse Gleichgewichtung.")
    print()
    print("Der Grund ist statistisch, nicht finanzwirtschaftlich: Um eine")
    print(f"Rendite von {WAHRE_RENDITE*100:.2f} % taeglich von null zu unterscheiden,")
    print(f"braucht man bei {WAHRE_SCHWANKUNG*100:.1f} % Schwankung ueber tausend")
    print("Beobachtungen. Mit 250 misst man ueberwiegend Rauschen - und der")
    print("Optimierer nimmt jedes Rauschen fuer bare Muenze.")
    print()
    print("Drei praktische Konsequenzen:")
    print("  1. Verzichten Sie auf Renditeschaetzungen, wo es geht. Das")
    print("     Minimum-Varianz-Portfolio braucht gar keine.")
    print("  2. Wenn Sie welche brauchen: schrumpfen Sie sie stark zur Mitte")
    print("     (James-Stein, Black-Litterman) - viel staerker als Kovarianzen.")
    print("  3. Begrenzen Sie Einzelgewichte. Was der Optimierer nicht darf,")
    print("     kann er auch nicht auf Rauschen setzen.")
    print("=" * 78)