{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Kapitel 13: Dynamische Programmierung — Die Bellman-Gleichung{idx:Bellman-Gleichung} und Order-Execution\n", "\n", "Begleitnotebook zu *Optimierte Entscheidungsfindung mit Python*. Die Codezellen sind identisch mit den im Buch abgedruckten Programmen.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Einmalig ausfuehren: installiert alle im Buch verwendeten Pakete.\n", "# Lokal in einer virtuellen Umgebung genauso gueltig wie in Google Colab.\n", "%pip install --quiet ortools highspy cvxpy scipy numpy pandas polars \\\n", " scikit-learn matplotlib plotly pyomo linopy pymoo pydantic openpyxl" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Die Bellman-Gleichung\n", "\n", "`Bellman_Minimalbeispiel.py`\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#!/usr/bin/env python3\n", "\n", "# Bellman_Minimalbeispiel.py\n", "\"\"\"\n", "Kapitel Dynamische Programmierung: Die Handrechnung zur Rueckwaertsinduktion als Code.\n", "Zeigt die Wertfunktionstabelle und die optimale Politik Schritt fuer Schritt.\n", "\"\"\"\n", "\n", "import numpy as np\n", "\n", "GESAMT = 3 # zu verkaufende Einheiten\n", "PERIODEN = 2 # Anzahl Verkaufsperioden\n", "\n", "\n", "def kosten(menge: int) -> float:\n", " \"\"\"Ueberproportionale Marktauswirkung: doppelte Menge kostet vierfach.\"\"\"\n", " return float(menge ** 2)\n", "\n", "\n", "def loese_rueckwaerts():\n", " # V[t, x] = minimale Restkosten, wenn zu Beginn von Periode t noch x Stueck offen sind\n", " V = np.full((PERIODEN + 1, GESAMT + 1), np.inf)\n", " politik = np.zeros((PERIODEN, GESAMT + 1), dtype=int)\n", "\n", " # Endbedingung: nach der letzten Periode darf nichts mehr offen sein\n", " V[PERIODEN, 0] = 0.0\n", "\n", " print(\"=\" * 70)\n", " print(\" RUECKWAERTSINDUKTION SCHRITT FUER SCHRITT\")\n", " print(\"=\" * 70)\n", "\n", " for t in range(PERIODEN - 1, -1, -1):\n", " letzte_periode = (t == PERIODEN - 1)\n", " print(f\"\\nStufe t = {t}\" + (\" (letzte Periode: alles muss weg)\" if letzte_periode\n", " else \" (freie Wahl der Menge)\"))\n", " print(f\" {'Zustand x':>10} | {'beste Aktion':>12} | {'Sofortkosten':>13} | \"\n", " f\"{'V[t+1]':>9} | {'V[t]':>8}\")\n", " print(\" \" + \"-\" * 62)\n", "\n", " for x in range(GESAMT + 1):\n", " aktionen = [x] if letzte_periode else range(x + 1)\n", " bester_wert, beste_aktion, beste_teile = np.inf, 0, (0.0, 0.0)\n", "\n", " for n in aktionen:\n", " rest = x - n\n", " sofort = kosten(n)\n", " zukunft = V[t + 1, rest]\n", " gesamt = sofort + zukunft\n", " if gesamt < bester_wert:\n", " bester_wert, beste_aktion = gesamt, n\n", " beste_teile = (sofort, zukunft)\n", "\n", " V[t, x] = bester_wert\n", " politik[t, x] = beste_aktion\n", " print(f\" {x:>10} | {beste_aktion:>12} | {beste_teile[0]:>13.1f} | \"\n", " f\"{beste_teile[1]:>9.1f} | {bester_wert:>8.1f}\")\n", "\n", " return V, politik\n", "\n", "\n", "if __name__ == \"__main__\":\n", " V, politik = loese_rueckwaerts()\n", "\n", " # --- Vorwaertspfad: der optimalen Politik folgen ---------------------\n", " print(\"\\n\" + \"=\" * 70)\n", " print(\" OPTIMALER PFAD (Vorwaertssimulation)\")\n", " print(\"=\" * 70)\n", " bestand = GESAMT\n", " gesamtkosten = 0.0\n", " for t in range(PERIODEN):\n", " aktion = politik[t, bestand]\n", " gesamtkosten += kosten(aktion)\n", " print(f\" Periode {t}: Bestand {bestand} -> verkaufe {aktion} \"\n", " f\"(Kosten {kosten(aktion):.1f}) -> Rest {bestand - aktion}\")\n", " bestand -= aktion\n", "\n", " print(f\"\\n Gesamtkosten: {gesamtkosten:.1f} (V[0, {GESAMT}] = {V[0, GESAMT]:.1f})\")\n", " assert abs(gesamtkosten - V[0, GESAMT]) < 1e-9, \"Pfadkosten != Wertfunktion!\"\n", "\n", " # --- Vergleich mit naiven Strategien ---------------------------------\n", " alles_sofort = kosten(GESAMT)\n", " print(f\"\\n Zum Vergleich - alles in Periode 0 verkaufen: {alles_sofort:.1f}\")\n", " print(f\" Ersparnis durch Stueckelung: {alles_sofort - gesamtkosten:.1f} \"\n", " f\"({(1 - gesamtkosten/alles_sofort)*100:.0f} %)\")\n", " print(\"=\" * 70)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Das Almgren-Chriss-Problem: optimale Orderausführung\n", "\n", "`Mehrperiodige_Order_Execution.py`\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#!/usr/bin/env python3\n", "\n", "# Mehrperiodige_Order_Execution.py\n", "\"\"\"\n", "Kapitel Dynamische Programmierung: Dynamische Programmierung fuer optimale Orderausfuehrung\n", "(Almgren-Chriss-Rahmen, geloest per Rueckwaertsinduktion).\n", "\n", "Eigenschaften:\n", " * Risikoterm sauber hergeleitet ueber den Aktienkurs P0 (Einheiten: EUR),\n", " ohne undokumentierte Skalierungsfaktoren\n", " * Vergleich mit der analytischen Almgren-Chriss-Loesung\n", " * Vergleich mit naiven Strategien (alles sofort / gleichmaessig)\n", " * Sensitivitaet gegenueber der Risikoaversion\n", "\"\"\"\n", "\n", "import os\n", "\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib\n", "matplotlib.use(\"Agg\")\n", "import matplotlib.pyplot as plt\n", "\n", "OUTPUT_DIR = os.path.join(os.path.dirname(os.path.abspath(__file__)), \"output\")\n", "os.makedirs(OUTPUT_DIR, exist_ok=True)\n", "\n", "# --- Parameter -------------------------------------------------------------\n", "GESAMTBESTAND = 100_000 # X_0, zu verkaufende Aktien\n", "PERIODEN = 5 # T Handelsperioden\n", "KURS = 50.0 # P_0 in EUR, zur Skalierung des Risikoterms\n", "ETA = 2.5e-6 # EUR je Stueck^2 (Slippage-Koeffizient)\n", "VOLA_JAHR = 0.30 # 30 % p.a.\n", "HANDELSSTUNDEN_JAHR = 252 * 6.5\n", "RISIKOAVERSION = 1e-5 # 1/EUR; kalibriert, siehe Kommentar unten\n", "SCHRITTWEITE = 1000 # Diskretisierung des Zustandsraums\n", "\n", "# Kalibrierungshinweis: Die dimensionslose Kennzahl des Modells ist\n", "# kappa_tilde^2 = lambda * sigma_periode^2 * P0^2 / eta\n", "# Sie entscheidet ueber den Charakter der Loesung:\n", "# << 1 -> praktisch gleichmaessige Aufteilung (Risiko spielt keine Rolle)\n", "# ~ 1 -> ausgewogener Kompromiss <- hier: 0.55\n", "# >> 1 -> fast alles sofort verkaufen\n", "# Genau diese Interpretierbarkeit geht mit einem undokumentierten\n", "# Skalierungsfaktor verloren.\n", "\n", "VOLA_PERIODE = VOLA_JAHR / np.sqrt(HANDELSSTUNDEN_JAHR)\n", "\n", "\n", "def periodenkosten(verkauf: float, restbestand: float) -> float:\n", " \"\"\"\n", " Sofortkosten einer Periode in EUR:\n", " (1) Marktauswirkung: eta * n^2\n", " (2) Risiko des Restbestands: lambda/2 * sigma^2 * P0^2 * X^2\n", " -> P0^2 macht aus \"Stueck^2\" einen EUR^2-Wert; lambda hat damit\n", " die Einheit 1/EUR und ist interpretierbar.\n", " \"\"\"\n", " marktauswirkung = ETA * verkauf ** 2\n", " wertvarianz = (VOLA_PERIODE ** 2) * (KURS ** 2) * (restbestand ** 2)\n", " risiko = 0.5 * RISIKOAVERSION * wertvarianz\n", " return marktauswirkung + risiko\n", "\n", "\n", "def loese_dp():\n", " \"\"\"Rueckwaertsinduktion ueber den diskretisierten Zustandsraum.\"\"\"\n", " zustaende = np.arange(0, GESAMTBESTAND + SCHRITTWEITE, SCHRITTWEITE)\n", " anzahl = len(zustaende)\n", "\n", " V = np.full((PERIODEN + 1, anzahl), np.inf)\n", " politik = np.zeros((PERIODEN, anzahl), dtype=int)\n", " V[PERIODEN, 0] = 0.0 # am Ende muss alles verkauft sein\n", "\n", " for t in range(PERIODEN - 1, -1, -1):\n", " for idx, bestand in enumerate(zustaende):\n", " if t == PERIODEN - 1:\n", " moegliche = [bestand] # letzte Periode: Rest muss weg\n", " else:\n", " moegliche = zustaende[zustaende <= bestand]\n", "\n", " bester_wert, beste_aktion = np.inf, 0\n", " for verkauf in moegliche:\n", " rest = bestand - verkauf\n", " rest_idx = int(round(rest / SCHRITTWEITE))\n", " gesamt = periodenkosten(verkauf, rest) + V[t + 1, rest_idx]\n", " if gesamt < bester_wert:\n", " bester_wert, beste_aktion = gesamt, int(verkauf)\n", "\n", " V[t, idx] = bester_wert\n", " politik[t, idx] = beste_aktion\n", "\n", " return zustaende, V, politik\n", "\n", "\n", "def analytische_loesung():\n", " \"\"\"\n", " Geschlossene Almgren-Chriss-Loesung fuer den kontinuierlichen Fall.\n", " Der optimale Pfad ist X_t = X_0 * sinh(kappa*(T-t)) / sinh(kappa*T)\n", " mit kappa = arccosh(tilde_kappa^2/2 + 1), tilde_kappa^2 = lambda*sigma^2*P0^2/eta.\n", " Dient hier als unabhaengige Kontrolle des DP-Ergebnisses.\n", " \"\"\"\n", " kappa_tilde_quadrat = (RISIKOAVERSION * (VOLA_PERIODE ** 2) * (KURS ** 2)) / ETA\n", " kappa = np.arccosh(kappa_tilde_quadrat / 2.0 + 1.0)\n", " if kappa < 1e-12: # Grenzfall: risikoneutral -> linear\n", " return np.linspace(GESAMTBESTAND, 0, PERIODEN + 1)\n", " t = np.arange(PERIODEN + 1)\n", " return GESAMTBESTAND * np.sinh(kappa * (PERIODEN - t)) / np.sinh(kappa * PERIODEN)\n", "\n", "\n", "def bewerte_pfad(bestaende):\n", " \"\"\"Gesamtkosten eines beliebigen Bestandspfades.\"\"\"\n", " summe = 0.0\n", " for t in range(len(bestaende) - 1):\n", " verkauf = bestaende[t] - bestaende[t + 1]\n", " summe += periodenkosten(verkauf, bestaende[t + 1])\n", " return summe\n", "\n", "\n", "if __name__ == \"__main__\":\n", " zustaende, V, politik = loese_dp()\n", "\n", " # --- Vorwaertspfad der optimalen Politik -----------------------------\n", " bestand = GESAMTBESTAND\n", " verlauf = [bestand]\n", " verkaeufe = []\n", " for t in range(PERIODEN):\n", " idx = int(round(bestand / SCHRITTWEITE))\n", " verkauf = politik[t, idx]\n", " verkaeufe.append(verkauf)\n", " bestand -= verkauf\n", " verlauf.append(bestand)\n", "\n", " print(\"=\" * 84)\n", " print(\" OPTIMALE MEHRPERIODIGE ORDER-EXECUTION (BELLMAN DP)\")\n", " print(\"=\" * 84)\n", " print(f\"Gesamtvolumen: {GESAMTBESTAND:,} Stueck zu {KURS:.2f} EUR \"\n", " f\"= {GESAMTBESTAND*KURS:,.0f} EUR Positionswert\")\n", " print(f\"Zeithorizont: {PERIODEN} Handelsperioden\")\n", " print(f\"Volatilitaet: {VOLA_JAHR*100:.0f} % p.a. \"\n", " f\"= {VOLA_PERIODE*100:.3f} % je Periode\")\n", " print(f\"Slippage eta: {ETA:.2e} EUR/Stueck^2\")\n", " print(f\"Risikoaversion: {RISIKOAVERSION:.2e} 1/EUR\")\n", " print(f\"Erwartete Gesamtreibung: {V[0, -1]:,.2f} EUR \"\n", " f\"({V[0, -1]/(GESAMTBESTAND*KURS)*10000:.1f} Basispunkte)\\n\")\n", "\n", " plan = pd.DataFrame([{\n", " \"Periode\": f\"t = {t} -> {t+1}\",\n", " \"Startbestand\": f\"{verlauf[t]:,}\",\n", " \"Verkauf n_t\": f\"{verkaeufe[t]:,}\",\n", " \"Restbestand\": f\"{verlauf[t+1]:,}\",\n", " \"Anteil\": f\"{verkaeufe[t]/GESAMTBESTAND*100:5.1f} %\",\n", " \"Kosten (EUR)\": f\"{periodenkosten(verkaeufe[t], verlauf[t+1]):,.0f}\",\n", " } for t in range(PERIODEN)])\n", " print(plan.to_string(index=False))\n", "\n", " # --- Vergleich mit Alternativen und der analytischen Loesung ---------\n", " sofort = [GESAMTBESTAND] + [0] * PERIODEN\n", " gleichmaessig = [GESAMTBESTAND * (1 - t / PERIODEN) for t in range(PERIODEN + 1)]\n", " analytisch = analytische_loesung()\n", "\n", " print(\"\\n\" + \"-\" * 84)\n", " print(f\"{'Strategie':<34} {'Kosten (EUR)':>15} {'Basispunkte':>13} \"\n", " f\"{'ggue. Optimum':>16}\")\n", " print(\"-\" * 84)\n", " optimum = V[0, -1]\n", " for name, pfad in [(\"DP-Optimum\", verlauf),\n", " (\"Analytisch (Almgren-Chriss)\", list(analytisch)),\n", " (\"Gleichmaessig (TWAP)\", gleichmaessig),\n", " (\"Alles sofort\", sofort)]:\n", " kosten = bewerte_pfad(pfad)\n", " bp = kosten / (GESAMTBESTAND * KURS) * 10000\n", " print(f\"{name:<34} {kosten:>15,.0f} {bp:>12.1f} \"\n", " f\"{kosten - optimum:>+15,.0f}\")\n", "\n", " print(\"-\" * 84)\n", " abweichung = abs(bewerte_pfad(list(analytisch)) - optimum) / optimum\n", " print(f\"Abweichung DP zur analytischen Loesung: {abweichung*100:.3f} % \"\n", " f\"(Diskretisierung: {SCHRITTWEITE} Stueck)\")\n", "\n", " # --- Sensitivitaet gegenueber der Risikoaversion ---------------------\n", " print(\"\\n--- Wie wirkt die Risikoaversion? ---\")\n", " print(f\"{'lambda':>10} | {'Verkauf in Periode 0':>22} | {'Charakter':<28}\")\n", " print(\"-\" * 70)\n", " for lam in [1e-7, 1e-6, 1e-5, 1e-4, 1e-3]:\n", " globals()[\"RISIKOAVERSION\"] = lam\n", " _, V_l, pol_l = loese_dp()\n", " erste = pol_l[0, -1]\n", " anteil = erste / GESAMTBESTAND * 100\n", " charakter = (\"nahezu gleichmaessig\" if anteil < 25 else\n", " \"front-loaded\" if anteil < 60 else \"fast alles sofort\")\n", " print(f\"{lam:>10.0e} | {erste:>13,} ({anteil:5.1f} %) | {charakter:<28}\")\n", " globals()[\"RISIKOAVERSION\"] = 1e-5 # zuruecksetzen\n", "\n", " # --- Diagramm ---------------------------------------------------------\n", " plt.figure(figsize=(10, 5.5))\n", " plt.plot(range(PERIODEN + 1), verlauf, \"o-\", lw=2.5, label=\"DP-Optimum\")\n", " plt.plot(range(PERIODEN + 1), analytisch, \"s--\", lw=1.8, alpha=0.8,\n", " label=\"Analytisch (Almgren-Chriss)\")\n", " plt.plot(range(PERIODEN + 1), gleichmaessig, \"^:\", lw=1.8, alpha=0.8,\n", " label=\"Gleichmaessig (TWAP)\")\n", " plt.bar(range(PERIODEN), verkaeufe, alpha=0.25, color=\"orange\", width=0.45,\n", " label=\"Verkaufstranche $n_t$\")\n", " plt.title(\"Optimaler Liquidationspfad ueber diskrete Perioden\", fontsize=12)\n", " plt.xlabel(\"Handelsperiode $t$\")\n", " plt.ylabel(\"Verbleibender Bestand $X_t$\")\n", " plt.xticks(range(PERIODEN + 1))\n", " plt.grid(True, linestyle=\":\", alpha=0.6)\n", " plt.legend()\n", " plt.tight_layout()\n", " ziel = os.path.join(OUTPUT_DIR, \"optimal_execution_dp.png\")\n", " plt.savefig(ziel, dpi=150)\n", " print(f\"\\nDiagramm gespeichert unter '{ziel}'\")\n", " print(\"=\" * 84)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.11" } }, "nbformat": 4, "nbformat_minor": 5 }