#!/usr/bin/env python3 # Mehrperiodige_Order_Execution.py """ Kapitel Dynamische Programmierung: Dynamische Programmierung fuer optimale Orderausfuehrung (Almgren-Chriss-Rahmen, geloest per Rueckwaertsinduktion). Eigenschaften: * Risikoterm sauber hergeleitet ueber den Aktienkurs P0 (Einheiten: EUR), ohne undokumentierte Skalierungsfaktoren * Vergleich mit der analytischen Almgren-Chriss-Loesung * Vergleich mit naiven Strategien (alles sofort / gleichmaessig) * Sensitivitaet gegenueber der Risikoaversion """ import os import numpy as np import pandas as pd import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt OUTPUT_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), "output") os.makedirs(OUTPUT_DIR, exist_ok=True) # --- Parameter ------------------------------------------------------------- GESAMTBESTAND = 100_000 # X_0, zu verkaufende Aktien PERIODEN = 5 # T Handelsperioden KURS = 50.0 # P_0 in EUR, zur Skalierung des Risikoterms ETA = 2.5e-6 # EUR je Stueck^2 (Slippage-Koeffizient) VOLA_JAHR = 0.30 # 30 % p.a. HANDELSSTUNDEN_JAHR = 252 * 6.5 RISIKOAVERSION = 1e-5 # 1/EUR; kalibriert, siehe Kommentar unten SCHRITTWEITE = 1000 # Diskretisierung des Zustandsraums # Kalibrierungshinweis: Die dimensionslose Kennzahl des Modells ist # kappa_tilde^2 = lambda * sigma_periode^2 * P0^2 / eta # Sie entscheidet ueber den Charakter der Loesung: # << 1 -> praktisch gleichmaessige Aufteilung (Risiko spielt keine Rolle) # ~ 1 -> ausgewogener Kompromiss <- hier: 0.55 # >> 1 -> fast alles sofort verkaufen # Genau diese Interpretierbarkeit geht mit einem undokumentierten # Skalierungsfaktor verloren. VOLA_PERIODE = VOLA_JAHR / np.sqrt(HANDELSSTUNDEN_JAHR) def periodenkosten(verkauf: float, restbestand: float) -> float: """ Sofortkosten einer Periode in EUR: (1) Marktauswirkung: eta * n^2 (2) Risiko des Restbestands: lambda/2 * sigma^2 * P0^2 * X^2 -> P0^2 macht aus "Stueck^2" einen EUR^2-Wert; lambda hat damit die Einheit 1/EUR und ist interpretierbar. """ marktauswirkung = ETA * verkauf ** 2 wertvarianz = (VOLA_PERIODE ** 2) * (KURS ** 2) * (restbestand ** 2) risiko = 0.5 * RISIKOAVERSION * wertvarianz return marktauswirkung + risiko def loese_dp(): """Rueckwaertsinduktion ueber den diskretisierten Zustandsraum.""" zustaende = np.arange(0, GESAMTBESTAND + SCHRITTWEITE, SCHRITTWEITE) anzahl = len(zustaende) V = np.full((PERIODEN + 1, anzahl), np.inf) politik = np.zeros((PERIODEN, anzahl), dtype=int) V[PERIODEN, 0] = 0.0 # am Ende muss alles verkauft sein for t in range(PERIODEN - 1, -1, -1): for idx, bestand in enumerate(zustaende): if t == PERIODEN - 1: moegliche = [bestand] # letzte Periode: Rest muss weg else: moegliche = zustaende[zustaende <= bestand] bester_wert, beste_aktion = np.inf, 0 for verkauf in moegliche: rest = bestand - verkauf rest_idx = int(round(rest / SCHRITTWEITE)) gesamt = periodenkosten(verkauf, rest) + V[t + 1, rest_idx] if gesamt < bester_wert: bester_wert, beste_aktion = gesamt, int(verkauf) V[t, idx] = bester_wert politik[t, idx] = beste_aktion return zustaende, V, politik def analytische_loesung(): """ Geschlossene Almgren-Chriss-Loesung fuer den kontinuierlichen Fall. Der optimale Pfad ist X_t = X_0 * sinh(kappa*(T-t)) / sinh(kappa*T) mit kappa = arccosh(tilde_kappa^2/2 + 1), tilde_kappa^2 = lambda*sigma^2*P0^2/eta. Dient hier als unabhaengige Kontrolle des DP-Ergebnisses. """ kappa_tilde_quadrat = (RISIKOAVERSION * (VOLA_PERIODE ** 2) * (KURS ** 2)) / ETA kappa = np.arccosh(kappa_tilde_quadrat / 2.0 + 1.0) if kappa < 1e-12: # Grenzfall: risikoneutral -> linear return np.linspace(GESAMTBESTAND, 0, PERIODEN + 1) t = np.arange(PERIODEN + 1) return GESAMTBESTAND * np.sinh(kappa * (PERIODEN - t)) / np.sinh(kappa * PERIODEN) def bewerte_pfad(bestaende): """Gesamtkosten eines beliebigen Bestandspfades.""" summe = 0.0 for t in range(len(bestaende) - 1): verkauf = bestaende[t] - bestaende[t + 1] summe += periodenkosten(verkauf, bestaende[t + 1]) return summe if __name__ == "__main__": zustaende, V, politik = loese_dp() # --- Vorwaertspfad der optimalen Politik ----------------------------- bestand = GESAMTBESTAND verlauf = [bestand] verkaeufe = [] for t in range(PERIODEN): idx = int(round(bestand / SCHRITTWEITE)) verkauf = politik[t, idx] verkaeufe.append(verkauf) bestand -= verkauf verlauf.append(bestand) print("=" * 84) print(" OPTIMALE MEHRPERIODIGE ORDER-EXECUTION (BELLMAN DP)") print("=" * 84) print(f"Gesamtvolumen: {GESAMTBESTAND:,} Stueck zu {KURS:.2f} EUR " f"= {GESAMTBESTAND*KURS:,.0f} EUR Positionswert") print(f"Zeithorizont: {PERIODEN} Handelsperioden") print(f"Volatilitaet: {VOLA_JAHR*100:.0f} % p.a. " f"= {VOLA_PERIODE*100:.3f} % je Periode") print(f"Slippage eta: {ETA:.2e} EUR/Stueck^2") print(f"Risikoaversion: {RISIKOAVERSION:.2e} 1/EUR") print(f"Erwartete Gesamtreibung: {V[0, -1]:,.2f} EUR " f"({V[0, -1]/(GESAMTBESTAND*KURS)*10000:.1f} Basispunkte)\n") plan = pd.DataFrame([{ "Periode": f"t = {t} -> {t+1}", "Startbestand": f"{verlauf[t]:,}", "Verkauf n_t": f"{verkaeufe[t]:,}", "Restbestand": f"{verlauf[t+1]:,}", "Anteil": f"{verkaeufe[t]/GESAMTBESTAND*100:5.1f} %", "Kosten (EUR)": f"{periodenkosten(verkaeufe[t], verlauf[t+1]):,.0f}", } for t in range(PERIODEN)]) print(plan.to_string(index=False)) # --- Vergleich mit Alternativen und der analytischen Loesung --------- sofort = [GESAMTBESTAND] + [0] * PERIODEN gleichmaessig = [GESAMTBESTAND * (1 - t / PERIODEN) for t in range(PERIODEN + 1)] analytisch = analytische_loesung() print("\n" + "-" * 84) print(f"{'Strategie':<34} {'Kosten (EUR)':>15} {'Basispunkte':>13} " f"{'ggue. Optimum':>16}") print("-" * 84) optimum = V[0, -1] for name, pfad in [("DP-Optimum", verlauf), ("Analytisch (Almgren-Chriss)", list(analytisch)), ("Gleichmaessig (TWAP)", gleichmaessig), ("Alles sofort", sofort)]: kosten = bewerte_pfad(pfad) bp = kosten / (GESAMTBESTAND * KURS) * 10000 print(f"{name:<34} {kosten:>15,.0f} {bp:>12.1f} " f"{kosten - optimum:>+15,.0f}") print("-" * 84) abweichung = abs(bewerte_pfad(list(analytisch)) - optimum) / optimum print(f"Abweichung DP zur analytischen Loesung: {abweichung*100:.3f} % " f"(Diskretisierung: {SCHRITTWEITE} Stueck)") # --- Sensitivitaet gegenueber der Risikoaversion --------------------- print("\n--- Wie wirkt die Risikoaversion? ---") print(f"{'lambda':>10} | {'Verkauf in Periode 0':>22} | {'Charakter':<28}") print("-" * 70) for lam in [1e-7, 1e-6, 1e-5, 1e-4, 1e-3]: globals()["RISIKOAVERSION"] = lam _, V_l, pol_l = loese_dp() erste = pol_l[0, -1] anteil = erste / GESAMTBESTAND * 100 charakter = ("nahezu gleichmaessig" if anteil < 25 else "front-loaded" if anteil < 60 else "fast alles sofort") print(f"{lam:>10.0e} | {erste:>13,} ({anteil:5.1f} %) | {charakter:<28}") globals()["RISIKOAVERSION"] = 1e-5 # zuruecksetzen # --- Diagramm --------------------------------------------------------- plt.figure(figsize=(10, 5.5)) plt.plot(range(PERIODEN + 1), verlauf, "o-", lw=2.5, label="DP-Optimum") plt.plot(range(PERIODEN + 1), analytisch, "s--", lw=1.8, alpha=0.8, label="Analytisch (Almgren-Chriss)") plt.plot(range(PERIODEN + 1), gleichmaessig, "^:", lw=1.8, alpha=0.8, label="Gleichmaessig (TWAP)") plt.bar(range(PERIODEN), verkaeufe, alpha=0.25, color="orange", width=0.45, label="Verkaufstranche $n_t$") plt.title("Optimaler Liquidationspfad ueber diskrete Perioden", fontsize=12) plt.xlabel("Handelsperiode $t$") plt.ylabel("Verbleibender Bestand $X_t$") plt.xticks(range(PERIODEN + 1)) plt.grid(True, linestyle=":", alpha=0.6) plt.legend() plt.tight_layout() ziel = os.path.join(OUTPUT_DIR, "optimal_execution_dp.png") plt.savefig(ziel, dpi=150) print(f"\nDiagramm gespeichert unter '{ziel}'") print("=" * 84)