operations_research/OR_HTML_04/Notebooks_04/dynamische-programmierung.ipynb

{
 "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": {
    "tags": [
     "nbmake-skip"
    ]
   },
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Notebook-spezifisches Setup: __file__ definieren (Programme nutzen es\n",
    "# fuer OUTPUT_DIR), Programmverzeichnis in den Suchpfad aufnehmen (or_kern)\n",
    "# und multiprocessing auf 'fork' stellen (Spawn findet Kernel-Funktionen nicht).\n",
    "import os, sys\n",
    "__file__ = os.path.join(os.getcwd(), '_notebook.py')\n",
    "for _p in [os.path.join(os.path.dirname(os.getcwd()),\n",
    "                       'Operations_Research_mit_Python_Version_04_Programme'),\n",
    "            os.path.dirname(os.getcwd()), os.getcwd()]:\n",
    "    if os.path.isfile(os.path.join(_p, 'or_kern.py')):\n",
    "        sys.path.insert(0, _p); break\n",
    "import multiprocessing as _mp\n",
    "try: _mp.set_start_method('fork', force=True)\n",
    "except (RuntimeError, ValueError): pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Die Bellman-Gleichung\n",
    "\n",
    "`Bellman_Minimalbeispiel.py`\n"
   ]
  },
  {
   "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": "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
}