#!/usr/bin/env python3 # CVaR_Portfolio.py """ Kapitel CVaR: CVaR-Optimierung mit L1-Transaktionskosten via CVXPY. Eigenschaften: * Spaltenreihenfolge erzwungen * Einheiten konsistent (alles taeglich, Annualisierung nur in der Ausgabe) * Szenario-Nebenbedingungen VEKTORISIERT statt in einer Python-Schleife (eine Matrixbedingung statt S einzelner Constraints - deutlich schneller) * Vergleich CVaR- gegen Varianz-Optimierung * Nachrechnung von VaR/CVaR aus den realisierten Szenarien """ import time import cvxpy as cp import numpy as np import pandas as pd import yfinance as yf TICKER = ["AAPL", "MSFT", "NVDA", "AMZN", "JNJ", "PFE", "JPM", "GS", "XOM", "CVX"] ALPHA = 0.95 # Konfidenzniveau: schlechteste 5 % der Tage MAX_GEWICHT = 0.25 GEBUEHRENSATZ = 0.002 # 0,2 % je Einheit Turnover (Spread + Brokerage) RISIKOAVERSION = 1.5 # bezogen auf TAEGLICHE Groessen HANDELSTAGE = 252 def lade_renditen(): ende = pd.Timestamp.today().normalize() start = ende - pd.DateOffset(years=2) roh = yf.download(TICKER, 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): kurse = roh["Close"][TICKER].dropna() # erzwingt eigene Spaltenreihenfolge else: kurse = roh[["Close"]].dropna() kurse.columns = TICKER assert list(kurse.columns) == TICKER, "Spaltenreihenfolge weicht ab!" return kurse.pct_change().dropna() def optimiere_cvar(R, w_alt, vektorisiert=True): """ Maximiere: taegliche Rendite - lambda * CVaR - Transaktionskosten Alle Groessen TAEGLICH. """ S, N = R.shape mu_taeglich = R.mean(axis=0) w = cp.Variable(N, nonneg=True) gamma = cp.Variable() # wird im Optimum zum VaR u = cp.Variable(S, nonneg=True) # Ueberschuss ueber die Schwelle cvar = gamma + (1.0 / (S * (1.0 - ALPHA))) * cp.sum(u) turnover = cp.norm1(w - w_alt) kosten = GEBUEHRENSATZ * turnover ziel = cp.Maximize(mu_taeglich @ w - RISIKOAVERSION * cvar - kosten) bedingungen = [cp.sum(w) == 1, w <= MAX_GEWICHT] if vektorisiert: # EINE Matrixbedingung statt S einzelner - deutlich schneller bedingungen.append(u >= -(R @ w) - gamma) else: for s in range(S): # Alternative: S einzelne Constraints (langsamer) bedingungen.append(u[s] >= -R[s] @ w - gamma) problem = cp.Problem(ziel, bedingungen) problem.solve() if problem.status not in ("optimal", "optimal_inaccurate"): raise SystemExit(f"CVaR-Problem nicht loesbar: {problem.status}") return w.value, float(gamma.value), float(cvar.value), float(turnover.value) def optimiere_varianz(R, w_alt): """Klassisches Mean-Variance zum Vergleich - ebenfalls taeglich gerechnet. ACHTUNG, DCP-Falle: Die Standardabweichung ist hier NICHT als cp.sqrt(cp.quad_form(w, sigma)) formulierbar. cp.sqrt ist konkav und verlangt ein konkaves Argument; quad_form ist konvex - CVXPY lehnt den Ausdruck mit einem DCPError ab, und zwar voellig zu Recht (Anhang Fehlerdiagnose). cp.psd_wrap() hilft dagegen nicht: Es behebt eine NUMERISCHE Beanstandung an sigma, keine Regelverletzung im Aufbau. Der Standardweg ist die Cholesky-Zerlegung sigma = L L^T. Damit gilt w' sigma w = ||L^T w||^2, also ist die Standardabweichung die 2-Norm eines AFFINEN Ausdrucks - konvex und damit regelkonform. """ N = R.shape[1] mu_taeglich = R.mean(axis=0) sigma = np.cov(R, rowvar=False, ddof=1) # Der winzige Diagonalzuschlag faengt den Fall ab, dass sigma numerisch # nur halbdefinit ist (mehr Titel als Handelstage, doppelte Spalten). L = np.linalg.cholesky(sigma + 1e-12 * np.eye(N)) w = cp.Variable(N, nonneg=True) ziel = cp.Maximize(mu_taeglich @ w - RISIKOAVERSION * cp.norm2(L.T @ w) - GEBUEHRENSATZ * cp.norm1(w - w_alt)) problem = cp.Problem(ziel, [cp.sum(w) == 1, w <= MAX_GEWICHT]) problem.solve() if problem.status not in ("optimal", "optimal_inaccurate"): raise SystemExit(f"Varianz-Problem nicht loesbar: {problem.status}") return w.value def realisierte_kennzahlen(w, R): """VaR und CVaR direkt aus den Szenarien - unabhaengige Gegenprobe.""" portfoliorenditen = R @ w verluste = -portfoliorenditen var = float(np.quantile(verluste, ALPHA)) cvar = float(verluste[verluste >= var].mean()) return var, cvar, float(portfoliorenditen.mean()), float(portfoliorenditen.std(ddof=1)) if __name__ == "__main__": renditen = lade_renditen() R = renditen.values S, N = R.shape w_alt = np.ones(N) / N # Ausgangslage: Gleichgewichtung t0 = time.perf_counter() w_cvar, var_modell, cvar_modell, turnover = optimiere_cvar(R, w_alt, vektorisiert=True) dauer_vektor = time.perf_counter() - t0 w_var = optimiere_varianz(R, w_alt) print("=" * 90) print(" CVaR-PORTFOLIO-OPTIMIERUNG MIT TRANSAKTIONSKOSTEN") print("=" * 90) print(f"Datenbasis: {S} Handelstage, {N} Titel | Konfidenzniveau " f"{ALPHA*100:.0f} % | Loesungszeit {dauer_vektor:.2f} s\n") # --- Gegenprobe: Modellwerte gegen realisierte Szenariowerte --------- var_real, cvar_real, mu_real, sd_real = realisierte_kennzahlen(w_cvar, R) print("--- Gegenprobe: stimmen Modell und Szenarien ueberein? ---") print(f" VaR aus dem Modell (gamma): {var_modell*100:7.4f} % | " f"aus den Szenarien: {var_real*100:7.4f} %") print(f" CVaR aus dem Modell: {cvar_modell*100:7.4f} % | " f"aus den Szenarien: {cvar_real*100:7.4f} %") assert abs(cvar_modell - cvar_real) < 1e-4, "CVaR stimmt nicht mit den Szenarien!" print(" -> Der Rockafellar-Uryasev-Trick liefert exakt den empirischen CVaR.") # --- Kennzahlen beider Portfolios ------------------------------------ print(f"\n{'Portfolio':<24} {'Rendite p.a.':>13} {'Vola p.a.':>11} " f"{'VaR 95% (Tag)':>15} {'CVaR 95% (Tag)':>16} {'Turnover':>10}") print("-" * 90) for name, w in [("CVaR-optimiert", w_cvar), ("Varianz-optimiert", w_var), ("Gleichgewichtung", w_alt)]: v, c, m, s = realisierte_kennzahlen(w, R) to = float(np.abs(w - w_alt).sum()) print(f"{name:<24} {m*HANDELSTAGE*100:>12.2f} % " f"{s*np.sqrt(HANDELSTAGE)*100:>10.2f} % {v*100:>14.3f} % " f"{c*100:>15.3f} % {to*100:>9.1f} %") print("\nHinweis zur Annualisierung: Renditen werden mit 252 skaliert,") print("Volatilitaeten mit sqrt(252). Fuer VaR/CVaR ist eine solche Skalierung") print("nur unter starken Annahmen (Unabhaengigkeit, kein Drift) zulaessig -") print("sie werden hier deshalb bewusst als TAGESwerte ausgewiesen.") # --- Allokationstabelle ------------------------------------------------ print("\n--- Allokation ---") tabelle = pd.DataFrame({ "Ticker": TICKER, "vorher": [f"{v*100:5.1f} %" for v in w_alt], "CVaR-opt.": [f"{v*100:5.1f} %" for v in w_cvar], "Handel": [f"{(w_cvar[i]-w_alt[i])*100:+6.1f} %" for i in range(N)], "Varianz-opt.": [f"{v*100:5.1f} %" for v in w_var], }) print(tabelle.to_string(index=False)) print(f"\nTurnover {turnover*100:.1f} % -> Transaktionskosten " f"{GEBUEHRENSATZ*turnover*100:.3f} % des Portfoliowerts") # --- Laufzeitvergleich vektorisiert vs. Schleife ---------------------- if S <= 600: # bei sehr vielen Szenarien zu langsam t0 = time.perf_counter() optimiere_cvar(R, w_alt, vektorisiert=False) dauer_schleife = time.perf_counter() - t0 print(f"\n--- Laufzeit: {S} Nebenbedingungen aufbauen ---") print(f" vektorisiert (u >= -(R @ w) - gamma): {dauer_vektor:6.2f} s") print(f" Schleife ueber Szenarien (V01): {dauer_schleife:6.2f} s " f"({dauer_schleife/dauer_vektor:.1f}x langsamer)") print("=" * 90)