operations_research/bilder_04/erzeuge_effizienzlinie.py

#!/usr/bin/env python3

# erzeuge_effizienzlinie.py
"""
Erzeugt die Figur "Was die Hausregeln kosten" zum Kapitel Markowitz:

    bilder_04/kap_markowitz_restriktionen.svg          statisch, fuer das PDF
    bilder_04/plotly/kap_markowitz_restriktionen.html  interaktiv, fuer die Website

Gezeigt werden zwei Effizienzlinien desselben Universums: einmal ohne
Beschraenkungen, einmal mit den institutionellen Regeln aus dem Kapitel
(hoechstens 20 % je Titel, Sektorgrenzen). Der senkrechte Abstand zwischen
beiden Kurven IST der Preis der Hausregeln - in Renditepunkten, bei gleichem
Risiko.

Warum eine feste Instanz und nicht die Kursdaten aus Markowitz_CVXPY.py:
Jenes Programm laedt Livedaten. Ein Diagramm daraus saehe an jedem Tag anders
aus, und der abgedruckte Text koennte nicht dazu passen. Hier wird deshalb ein
FESTES Universum aus einem Faktormodell erzeugt (Marktfaktor, vier
Sektorfaktoren, titelspezifisches Rauschen) - dieselbe Struktur wie im Kapitel,
dieselben Regeln, aber reproduzierbar. Die Zahlen des Livelaufs weichen ab, die
Aussage nicht.

Die interaktive Fassung zeigt im Hover, was das statische Bild nicht kann: die
GEWICHTE des jeweiligen Portfolios. Damit wird sichtbar, wie sich die
Zusammensetzung entlang der Linie verschiebt - und woran die beschraenkte
Fassung scheitert, sobald ein einzelner Titel an seine 20-%-Grenze stoesst.

Aufruf (aus dem Repository-Wurzelverzeichnis):
    python3 bilder_04/erzeuge_effizienzlinie.py

Benoetigt: numpy, cvxpy, matplotlib, plotly
"""

from __future__ import annotations

import os

import cvxpy as cp
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import plotly.graph_objects as go

# Reproduzierbare SVG-Ausgabe (siehe erzeuge_titelseite.py)
plt.rcParams["svg.hashsalt"] = "or-mit-python-v04"

BASIS = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
BILDER = os.path.join(BASIS, "bilder_04")
PLOTLY_VERZ = os.path.join(BILDER, "plotly")

INDIGO, CYAN, AMBER, GRUEN = "#4338ca", "#0891b2", "#b45309", "#059669"

# Struktur wie in Markowitz_CVXPY.py - zehn Titel in vier Sektoren
TITEL = ["AAPL", "MSFT", "NVDA", "AMZN", "JNJ", "PFE", "JPM", "GS", "XOM", "CVX"]
SEKTOR_VON = ["Technologie"] * 4 + ["Gesundheit"] * 2 + ["Finanzen"] * 2 + ["Energie"] * 2
SEKTORGRENZEN = {"Technologie": 0.35, "Gesundheit": 0.40,
                 "Finanzen": 0.40, "Energie": 0.40}
MAX_GEWICHT = 0.20
PUNKTE = 40


def universum(saat: int = 3) -> tuple[np.ndarray, np.ndarray]:
    """Erwartete Renditen und Kovarianzmatrix aus einem Faktormodell.

    Jeder Titel haengt am Markt, an seinem Sektor und an eigenem Rauschen.
    Genau diese Struktur macht Diversifikation ueberhaupt moeglich: Titel
    desselben Sektors laufen staerker im Gleichschritt als Titel verschiedener
    Sektoren - und nur der zweite Fall traegt.
    """
    rng = np.random.default_rng(saat)
    n = len(TITEL)
    sektoren = sorted(set(SEKTOR_VON))
    ladung_markt = rng.uniform(0.7, 1.3, n)
    ladung_sektor = rng.uniform(0.4, 0.9, n)
    eigen = rng.uniform(0.10, 0.26, n)

    sigma_markt, sigma_sektor = 0.16, 0.12
    kovarianz = np.outer(ladung_markt, ladung_markt) * sigma_markt ** 2
    for sektor in sektoren:
        maske = np.array([s == sektor for s in SEKTOR_VON], dtype=float)
        gewicht = ladung_sektor * maske
        kovarianz += np.outer(gewicht, gewicht) * sigma_sektor ** 2
    kovarianz += np.diag(eigen ** 2)

    # Erwartete Renditen: Marktprämie plus etwas titelspezifische Streuung
    mu = 0.03 + 0.055 * ladung_markt + rng.normal(0, 0.012, n)
    return mu, kovarianz


def sektormasken() -> dict[str, np.ndarray]:
    return {sektor: np.array([s == sektor for s in SEKTOR_VON], dtype=float)
            for sektor in SEKTORGRENZEN}


def linie(mu: np.ndarray, kovarianz: np.ndarray, beschraenkt: bool):
    """Effizienzlinie ueber das epsilon-Constraint-Verfahren.

    Fuer jede Zielrendite wird das minimale Risiko gesucht. Das ist dasselbe
    Verfahren, mit dem im Kapitel Mehrziel die Pareto-Front entsteht - eine
    Effizienzlinie IST eine Pareto-Front aus Rendite und Risiko.
    """
    n = len(mu)
    w = cp.Variable(n, nonneg=True)
    ziel_rendite = cp.Parameter(nonneg=True)

    # Cholesky statt quad_form: die 2-Norm eines affinen Ausdrucks ist
    # zweifelsfrei DCP (siehe die Fehlerdiagnose zu DCPError).
    wurzel = np.linalg.cholesky(kovarianz + 1e-12 * np.eye(n))

    bedingungen = [cp.sum(w) == 1, mu @ w >= ziel_rendite]
    if beschraenkt:
        bedingungen.append(w <= MAX_GEWICHT)
        for sektor, maske in sektormasken().items():
            bedingungen.append(maske @ w <= SEKTORGRENZEN[sektor])

    problem = cp.Problem(cp.Minimize(cp.norm2(wurzel.T @ w)), bedingungen)

    # Der obere Rand ist NICHT max(mu): Unter Restriktionen ist der beste
    # Einzeltitel gar nicht voll kaufbar. Also wird er ausgerechnet.
    w_max = cp.Variable(n, nonneg=True)
    obere = [cp.sum(w_max) == 1]
    if beschraenkt:
        obere.append(w_max <= MAX_GEWICHT)
        for sektor, maske in sektormasken().items():
            obere.append(maske @ w_max <= SEKTORGRENZEN[sektor])
    cp.Problem(cp.Maximize(mu @ w_max), obere).solve()
    hoechste = float(mu @ w_max.value)

    # Untere Grenze: die Rendite des Minimum-Varianz-Portfolios
    ziel_rendite.value = 0.0
    problem.solve()
    kleinste = float(mu @ w.value)

    punkte = []
    for ziel in np.linspace(kleinste, hoechste, PUNKTE):
        ziel_rendite.value = float(ziel)
        problem.solve()
        if problem.status not in (cp.OPTIMAL, cp.OPTIMAL_INACCURATE):
            continue
        gewichte = np.asarray(w.value).ravel()
        risiko = float(np.sqrt(gewichte @ kovarianz @ gewichte))
        punkte.append((risiko, float(mu @ gewichte), gewichte))
    return punkte


def aufpreis(frei, gebunden) -> tuple[float, float]:
    """Groesster senkrechter Abstand: Wo kosten die Regeln am meisten Rendite?

    Verglichen wird bei GLEICHEM Risiko - alles andere waere kein Preis,
    sondern nur ein anderer Punkt auf der Kurve.
    """
    risiken_frei = np.array([p[0] for p in frei])
    renditen_frei = np.array([p[1] for p in frei])
    bester_abstand, bestes_risiko = 0.0, 0.0
    for risiko, rendite, _ in gebunden:
        if risiko < risiken_frei.min() or risiko > risiken_frei.max():
            continue
        vergleich = float(np.interp(risiko, risiken_frei, renditen_frei))
        if vergleich - rendite > bester_abstand:
            bester_abstand, bestes_risiko = vergleich - rendite, risiko
    return bester_abstand, bestes_risiko


def zeichne_statisch(frei, gebunden, abstand, risiko_stelle) -> None:
    figur, achse = plt.subplots(figsize=(7.4, 4.6))
    for punkte, farbe, name, stil in (
            (frei, INDIGO, "ohne Beschränkungen", "-"),
            (gebunden, AMBER, "mit Hausregeln (20 % je Titel, Sektorgrenzen)", "-")):
        x = [p[0] * 100 for p in punkte]
        y = [p[1] * 100 for p in punkte]
        achse.plot(x, y, stil, color=farbe, linewidth=2.0, label=name)

    if abstand > 0:
        oben = float(np.interp(risiko_stelle,
                               [p[0] for p in frei], [p[1] for p in frei]))
        unten = oben - abstand
        achse.annotate("", xy=(risiko_stelle * 100, oben * 100),
                       xytext=(risiko_stelle * 100, unten * 100),
                       arrowprops=dict(arrowstyle="<->", color=GRUEN, lw=1.8))
        text = f"{abstand * 100:.2f} Prozentpunkte".replace(".", ",")
        achse.annotate(f"{text}\nRendite, bei gleichem Risiko",
                       (risiko_stelle * 100, (oben + unten) / 2 * 100),
                       textcoords="offset points", xytext=(10, 0),
                       fontsize=8.5, color=GRUEN, va="center")

    # Der teuerste Effekt ist nicht der senkrechte Abstand, sondern die Decke:
    # Die beschraenkte Linie endet frueher. Alles rechts davon ist mit den
    # Hausregeln ueberhaupt nicht erreichbar - egal, wie viel Risiko man traegt.
    ende_gebunden, ende_frei = gebunden[-1], frei[-1]
    achse.scatter([ende_gebunden[0] * 100], [ende_gebunden[1] * 100], s=110,
                  facecolors="none", edgecolors=AMBER, linewidths=2.2, zorder=5)
    verlust = (ende_frei[1] - ende_gebunden[1]) * 100
    achse.annotate(f"hier ist Schluss: {ende_gebunden[1] * 100:.2f} % statt "
                   f"{ende_frei[1] * 100:.2f} %\n"
                   f"— {verlust:.2f} Punkte, die kein Risiko mehr kauft".replace(".", ","),
                   (ende_gebunden[0] * 100, ende_gebunden[1] * 100),
                   textcoords="offset points", xytext=(8, -22),
                   fontsize=8.5, color=AMBER)

    achse.set_xlabel("Volatilität $\\sigma$ (% p. a.)")
    achse.set_ylabel("erwartete Rendite $\\mu$ (% p. a.)")
    achse.set_title("Was die Hausregeln kosten", fontsize=11)
    achse.grid(linestyle=":", alpha=0.5)
    achse.legend(fontsize=8.5, loc="lower right", frameon=False)
    for rand in ("top", "right"):
        achse.spines[rand].set_visible(False)
    figur.tight_layout()

    os.makedirs(BILDER, exist_ok=True)
    for endung in ("svg",):
        pfad = os.path.join(BILDER, f"kap_markowitz_restriktionen.{endung}")
        figur.savefig(pfad, format=endung, dpi=160,
                      metadata={"Date": None} if endung == "svg" else None)
        print(f"geschrieben: {pfad}")
    plt.close(figur)


def zeichne_interaktiv(frei, gebunden) -> None:
    figur = go.Figure()

    def beschriftung(punkt):
        risiko, rendite, gewichte = punkt
        reihenfolge = np.argsort(-gewichte)[:4]
        anteile = "<br>".join(
            f"&nbsp;&nbsp;{TITEL[i]} ({SEKTOR_VON[i]}): {gewichte[i] * 100:.1f} %"
            for i in reihenfolge if gewichte[i] > 0.005)
        sektoren = "<br>".join(
            f"&nbsp;&nbsp;{sektor}: {float(maske @ gewichte) * 100:.1f} %"
            for sektor, maske in sektormasken().items())
        return (f"<b>{rendite * 100:.2f} % Rendite</b><br>"
                f"{risiko * 100:.2f} % Volatilität<br>"
                f"Sharpe {(rendite - 0.03) / risiko:.2f}<br><br>"
                f"<b>größte Positionen</b><br>{anteile}<br><br>"
                f"<b>je Sektor</b><br>{sektoren}<extra></extra>")

    for punkte, farbe, name in ((frei, INDIGO, "ohne Beschränkungen"),
                                (gebunden, AMBER, "mit Hausregeln")):
        figur.add_trace(go.Scatter(
            x=[p[0] * 100 for p in punkte], y=[p[1] * 100 for p in punkte],
            mode="lines+markers", name=name,
            line=dict(color=farbe, width=2.2), marker=dict(size=6),
            hovertemplate=[beschriftung(p) for p in punkte]))

    figur.update_layout(
        title="Effizienzlinie mit und ohne Hausregeln — Gewichte im Hover",
        xaxis=dict(title="Volatilität σ (% p. a.)"),
        yaxis=dict(title="erwartete Rendite μ (% p. a.)"),
        height=460, margin=dict(l=70, r=20, t=55, b=50),
        paper_bgcolor="rgba(0,0,0,0)", plot_bgcolor="rgba(0,0,0,0)",
        legend=dict(orientation="h", y=-0.18), hovermode="closest")

    os.makedirs(PLOTLY_VERZ, exist_ok=True)
    pfad = os.path.join(PLOTLY_VERZ, "kap_markowitz_restriktionen.html")
    with open(pfad, "w", encoding="utf-8") as datei:
        datei.write(figur.to_html(full_html=False, include_plotlyjs=False,
                                  div_id="kap-markowitz-restriktionen"))
    print(f"geschrieben: {pfad}")


if __name__ == "__main__":
    mu, kovarianz = universum()
    frei = linie(mu, kovarianz, beschraenkt=False)
    gebunden = linie(mu, kovarianz, beschraenkt=True)
    abstand, stelle = aufpreis(frei, gebunden)
    print(f"frei: {len(frei)} Punkte, Rendite "
          f"{frei[0][1] * 100:.2f} bis {frei[-1][1] * 100:.2f} %")
    print(f"gebunden: {len(gebunden)} Punkte, Rendite "
          f"{gebunden[0][1] * 100:.2f} bis {gebunden[-1][1] * 100:.2f} %")
    print(f"groesster Abstand: {abstand * 100:.2f} Prozentpunkte "
          f"bei sigma = {stelle * 100:.2f} %")
    zeichne_statisch(frei, gebunden, abstand, stelle)
    zeichne_interaktiv(frei, gebunden)