#!/usr/bin/env python3

# VaR_CVaR_Demo.py
"""
Kapitel CVaR: Fat Tails, VaR und CVaR anschaulich.

Teil 1: Wie oft treten "unmoegliche" Tage wirklich auf?
Teil 2: Warum ist der VaR nicht subadditiv - ein Gegenbeispiel zum Nachrechnen.
"""

import numpy as np
from scipy import stats


def var_quantil(verluste: np.ndarray, alpha: float = 0.95) -> float:
    """VaR = Quantil der Verlustverteilung (Verluste positiv, Gewinne negativ)."""
    return float(np.quantile(verluste, alpha))


def cvar_rockafellar(verluste: np.ndarray, alpha: float = 0.95) -> float:
    """
    CVaR ueber die Rockafellar-Uryasev-Formel:
        CVaR = min_gamma { gamma + 1/(1-alpha) * E[max(Verlust - gamma, 0)] }

    WICHTIG: Der naheliegende Weg "Mittelwert aller Werte >= VaR" ist FALSCH,
    sobald die Verteilung Atome hat (z. B. genau zwei moegliche Verluste).
    Dann liegt der VaR selbst auf einem Atom, und der Vergleich '>=' erfasst
    zu viel Wahrscheinlichkeitsmasse. Die Formel unten behandelt das korrekt -
    und ist zugleich genau der Ausdruck, den wir im Abschnitt 'Value at Risk
    und Conditional Value at Risk' optimieren.
    """
    kandidaten = np.unique(verluste)          # Optimum liegt immer auf einem Datenpunkt
    return float(min(g + np.mean(np.maximum(verluste - g, 0.0)) / (1.0 - alpha)
                     for g in kandidaten))


if __name__ == "__main__":
    rng = np.random.default_rng(2026)

    # --- Teil 1: Fat Tails (analytisch, nicht simuliert) ------------------
    print("=" * 88)
    print("  TEIL 1: WIE OFT TRITT DAS 'UNMOEGLICHE' EIN?")
    print("=" * 88)
    print("Vergleich: Normalverteilung gegen t-Verteilung mit 3 Freiheitsgraden")
    print("(beide auf Standardabweichung 1 normiert).\n")

    t_verteilung = stats.t(df=3)
    skalierung = t_verteilung.std()            # auf Varianz 1 bringen

    print(f"{'Ereignis':<22} {'Normal':>14} {'t (df=3)':>14} {'Faktor':>11} "
          f"{'Normal: 1 Tag in':>18}")
    print("-" * 88)
    for k in [3, 4, 5, 6]:
        p_normal = 2 * stats.norm.sf(k)                     # beidseitig
        p_t = 2 * t_verteilung.sf(k * skalierung)
        print(f"Abweichung > {k} Sigma  {p_normal*100:>13.6f} % {p_t*100:>13.6f} % "
              f"{p_t/p_normal:>10.1f}x {1/p_normal/252:>15,.0f} Jahre")

    print("\nDeutung: Ein 5-Sigma-Tag ist unter Normalverteilung ein Ereignis von")
    print("etwa einmal in 6.900 Jahren. Reale Aktienmaerkte liefern solche Tage")
    print("mehrfach pro Jahrzehnt. Wer allein mit Varianz steuert, plant fuer")
    print("eine Welt, in der Crashs praktisch nicht vorkommen.")

    # --- Teil 2: VaR ist nicht subadditiv ---------------------------------
    print("\n" + "=" * 88)
    print("  TEIL 2: WARUM DER VaR KEIN KOHAERENTES RISIKOMASS IST")
    print("=" * 88)
    print("Zwei unabhaengige Anleihen, je 100 EUR Nominal.")
    print("Jede faellt mit 4 % Wahrscheinlichkeit aus (Verlust 100),")
    print("sonst zahlt sie 2 EUR Kupon (Verlust -2).\n")

    ziehungen = 2_000_000
    verlust_a = np.where(rng.random(ziehungen) < 0.04, 100.0, -2.0)
    verlust_b = np.where(rng.random(ziehungen) < 0.04, 100.0, -2.0)
    verlust_ab = verlust_a + verlust_b

    print(f"{'':<28} {'VaR 95%':>12} {'CVaR 95%':>12}")
    print("-" * 88)
    werte = {}
    for name, v in [("Anleihe A allein", verlust_a),
                    ("Anleihe B allein", verlust_b),
                    ("Portfolio A+B", verlust_ab)]:
        werte[name] = (var_quantil(v), cvar_rockafellar(v))
        print(f"{name:<28} {werte[name][0]:>12.2f} {werte[name][1]:>12.2f}")

    var_summe = werte["Anleihe A allein"][0] + werte["Anleihe B allein"][0]
    cvar_summe = werte["Anleihe A allein"][1] + werte["Anleihe B allein"][1]
    print(f"{'Summe der Einzelwerte':<28} {var_summe:>12.2f} {cvar_summe:>12.2f}")

    var_port, cvar_port = werte["Portfolio A+B"]
    print("-" * 88)
    print(f"VaR:  Portfolio {var_port:7.2f}  vs. Summe {var_summe:7.2f}  -> "
          f"{'VERLETZT die Subadditivitaet!' if var_port > var_summe else 'subadditiv'}")
    print(f"CVaR: Portfolio {cvar_port:7.2f}  vs. Summe {cvar_summe:7.2f}  -> "
          f"{'subadditiv (kohaerent)' if cvar_port <= cvar_summe + 1e-6 else 'verletzt'}")

    print("\nDeutung: Einzeln betrachtet meldet der VaR fuer jede Anleihe einen")
    print("GEWINN von 2 EUR - denn mit 96 % Wahrscheinlichkeit passiert nichts,")
    print("und 4 % liegen unterhalb der 5-%-Schwelle. Im Portfolio steigt die")
    print("Wahrscheinlichkeit mindestens eines Ausfalls auf 7,8 % und damit UEBER")
    print("die Schwelle - der VaR springt auf 98. Er behauptet also, Streuung")
    print("habe das Risiko erhoeht. Das ist oekonomisch unsinnig und der Grund,")
    print("warum die Bankenaufsicht mit Basel III auf den Expected Shortfall")
    print("umgestellt hat.")
    print("=" * 88)