#!/usr/bin/env python3 # Strukturbruecke.py """ Kapitel Bruecke: Derselbe Code, zwei Welten - der Beweis statt der Behauptung. Der Teil-Auftakt stellt eine Tabelle auf: "Ressourcen auf Produkte verteilen" entspreche "Kapital auf Anlagen verteilen", "gegen den Worst Case absichern" entspreche "Absicherung gegen Kursabstuerze". Solche Tabellen stehen in vielen Buechern. Sie sind billig - und man kann sie pruefen. Dieses Programm prueft sie. Es fuettert zweimal DENSELBEN Code mit Daten aus zwei Welten und zeigt die Ergebnisse nebeneinander: 1. Die Allokation als LP. Das Domaenenmodell aus or_kern.py, einmal mit einer Schreinerei (Montagestunden, Plattenmaterial) und einmal mit einem Depot (Kapital, Risikobudget). Gleiche Klasse, gleicher Modellbauer, gleiche Abnahmepruefung - und der Schattenpreis heisst in der einen Welt "Wert einer zusaetzlichen Montagestunde" und in der anderen "Preis des Risikos". 2. Die Absicherung gegen den schlechtesten Fall als CVaR. EINE Funktion, einmal mit Kursrenditen und einmal mit Lieferverzuegen. 3. Was NICHT hinueberreicht. Der ehrliche Teil: Drei Unterschiede, die die Analogie begrenzen - und die man kennen muss, bevor man Methoden aus der einen Welt in die andere traegt. ZUR SOLVERWAHL: Teil 1 benutzt 'loese_mit_scipy' und nicht 'loese_mit_glop'. Das ist kein Zufall - CVXPY laedt fuer Teil 2 highspy, und ortools vertraegt sich damit nicht im selben Prozess (Kapitel Oekosystem). Wer hier GLOP nimmt, bekommt beim cvxpy-Import eine Fehlermeldung ueber ein 'undefined symbol'. Genau deshalb laedt or_kern.py seine Solver erst beim Aufruf. Benoetigt: numpy, cvxpy, scipy und pydantic (ueber or_kern) """ from __future__ import annotations import numpy as np from or_kern import (Produkt, Produktionsproblem, loese_mit_scipy, pruefe_loesung) SAAT = 7 SZENARIEN = 500 ALPHA = 0.95 # die schlechtesten 5 % der Faelle MAX_ANTEIL = 0.40 # Streuungsgebot in beiden Welten # --- Teil 1: Dieselbe Klasse, zwei Welten ---------------------------------- def schreinerei() -> Produktionsproblem: """Montagestunden und Plattenmaterial auf Tische und Stuehle verteilen.""" return Produktionsproblem( produkte=[ Produkt(name="Tisch", deckungsbeitrag=240.0, verbrauch={"Montagestunden": 3.0, "Plattenmaterial": 6.0}), Produkt(name="Stuhl", deckungsbeitrag=60.0, verbrauch={"Montagestunden": 1.0, "Plattenmaterial": 1.0}), Produkt(name="Regal", deckungsbeitrag=130.0, verbrauch={"Montagestunden": 2.0, "Plattenmaterial": 4.0}), ], kapazitaeten={"Montagestunden": 150.0, "Plattenmaterial": 240.0}) def depot() -> Produktionsproblem: """Kapital und Risikobudget auf Anlageklassen verteilen. Eine "Einheit" ist hier 1.000 EUR Anlagesumme. Der Deckungsbeitrag ist der erwartete Jahresertrag dieser Einheit, der Verbrauch die beanspruchte Kapital- und Risikomenge. Das ist keine Analogie, sondern buchstaeblich dasselbe Modell - deshalb passt es in dieselbe Klasse. """ return Produktionsproblem( produkte=[ Produkt(name="Aktien Welt", deckungsbeitrag=75.0, verbrauch={"Kapital (Tsd. EUR)": 1.0, "Risikobudget": 1.00}), Produkt(name="Anleihen", deckungsbeitrag=28.0, verbrauch={"Kapital (Tsd. EUR)": 1.0, "Risikobudget": 0.22}), Produkt(name="Immobilienfonds", deckungsbeitrag=46.0, verbrauch={"Kapital (Tsd. EUR)": 1.0, "Risikobudget": 0.55}), ], kapazitaeten={"Kapital (Tsd. EUR)": 150.0, "Risikobudget": 90.0}) def berichte_allokation(titel: str, problem: Produktionsproblem, einheit: str, ertragsname: str) -> None: """Ein Bericht fuer beide Welten - nur die Beschriftung wechselt.""" loesung = loese_mit_scipy(problem) beanstandungen = pruefe_loesung(problem, loesung) print(f" {titel}") print(f" {loesung.als_bericht()}") for produkt in problem.produkte: print(f" {produkt.name:<18} {loesung.werte[produkt.name]:8.2f} {einheit}") print(f" {ertragsname:<18} {loesung.zielwert:8.2f} EUR") for ressource, preis in loesung.schattenpreise.items(): print(f" Schattenpreis {ressource:<24} {preis:7.2f} EUR") print(f" Abnahmepruefung: " f"{'bestanden' if not beanstandungen else beanstandungen}") # --- Teil 2: Eine CVaR-Funktion, zwei Welten ------------------------------- def optimiere_cvar(verluste: np.ndarray, ertrag: np.ndarray, mindestertrag: float, alpha: float = ALPHA, max_anteil: float = MAX_ANTEIL): """Minimiert den CVaR der Verluste unter einer Mindestertragsbedingung. 'verluste' hat die Form (Szenarien x Optionen) und enthaelt, was in Szenario s passiert, wenn eine Einheit in Option j steckt. Was ein "Verlust" ist, entscheidet allein die Einheit der Matrix: Prozentpunkte Kursverlust oder Tage Lieferverzug - die Formel sieht keinen Unterschied. Das ist die Rockafellar-Uryasev-Formulierung aus dem Kapitel CVaR, hier ohne jede Aenderung wiederverwendet. """ import cvxpy as cp anzahl_szenarien, anzahl_optionen = verluste.shape anteil = cp.Variable(anzahl_optionen, nonneg=True) schwelle = cp.Variable() # wird im Optimum zum VaR ueberschuss = cp.Variable(anzahl_szenarien, nonneg=True) cvar = schwelle + (1.0 / (anzahl_szenarien * (1 - alpha))) * cp.sum(ueberschuss) problem = cp.Problem( cp.Minimize(cvar), [ueberschuss >= verluste @ anteil - schwelle, cp.sum(anteil) == 1, anteil <= max_anteil, ertrag @ anteil >= mindestertrag]) problem.solve() if problem.status not in ("optimal", "optimal_inaccurate"): raise SystemExit(f"CVaR-Problem nicht loesbar: {problem.status}") return anteil.value, float(cvar.value), float(schwelle.value) def kursszenarien(rng) -> tuple[np.ndarray, np.ndarray, list[str]]: """Taegliche Verluste (negative Renditen) von sechs Anlageklassen.""" namen = ["Aktien Welt", "Aktien EU", "Schwellenlaender", "Staatsanleihen", "Unternehmensanl.", "Rohstoffe"] rendite_pa = np.array([0.080, 0.065, 0.110, 0.025, 0.045, 0.070]) schwankung = np.array([0.180, 0.160, 0.260, 0.040, 0.075, 0.210]) taeglich = rng.normal(rendite_pa / 252, schwankung / np.sqrt(252), (SZENARIEN, len(namen))) return -taeglich * 100.0, rendite_pa * 100.0, namen def lieferszenarien(rng) -> tuple[np.ndarray, np.ndarray, list[str]]: """Lieferverzug in Tagen bei sechs Lieferanten. Der Aufbau ist bewusst anders als bei den Kursen: Hier gibt es einen normalen Verzug UND seltene Totalausfaelle, die 14 Tage kosten. Das ist genau die Art fetter Raender, wegen der man in beiden Welten CVaR statt Standardabweichung benutzt. """ namen = ["Nordwerk", "Sued-Metall", "Fernost A", "Lokalzulieferer", "Fernost B", "Osteuropa"] zuverlaessigkeit = np.array([0.94, 0.97, 0.90, 0.99, 0.92, 0.95]) verzugsstreuung = np.array([3.5, 1.8, 5.0, 0.9, 4.2, 2.8]) normaler_verzug = np.maximum(0.0, rng.normal(0.5, 1.0, (SZENARIEN, len(namen))) * verzugsstreuung) totalausfall = (rng.random((SZENARIEN, len(namen))) > zuverlaessigkeit) * 14.0 return normaler_verzug + totalausfall, zuverlaessigkeit, namen def berichte_cvar(titel: str, verluste, ertrag, namen, mindestertrag, einheit: str, ertragsname: str): anteile, cvar, var = optimiere_cvar(verluste, ertrag, mindestertrag) mittel = float((verluste @ anteile).mean()) print(f" {titel}") print(f" VaR {ALPHA:.0%}: {var:8.3f} {einheit}") print(f" CVaR {ALPHA:.0%}: {cvar:8.3f} {einheit} " f"(Mittelwert ueber alle Szenarien: {mittel:.3f})") print(f" {ertragsname}: {float(ertrag @ anteile):.3f}") print(" Aufteilung:") for name, anteil in zip(namen, anteile): balken = "#" * int(round(anteil * 40)) grenze = " <- an der Streuungsgrenze" if anteil > MAX_ANTEIL - 1e-4 else "" print(f" {name:<18} {anteil:6.1%} {balken}{grenze}") return anteile if __name__ == "__main__": rng = np.random.default_rng(SAAT) print("=" * 80) print(" DIESELBE STRUKTUR, ZWEI WELTEN") print("=" * 80) # --- Teil 1 ---------------------------------------------------------- print("\n1. Allokation als LP - EINE Klasse, EIN Modellbauer\n") berichte_allokation("Werkstatt: Produktionsprogramm", schreinerei(), "Stueck", "Deckungsbeitrag") print() berichte_allokation("Depot: Anlageaufteilung", depot(), "Tsd. ", "Erwarteter Ertrag") print("\n Beide Ausgaben stammen aus derselben Funktion " "'berichte_allokation'.") print(" Ausgetauscht wurden nur die Daten und die Beschriftungen -") print(" keine Zeile Modellcode.") print() print(" Lesen Sie die Schattenpreise nebeneinander: In der Werkstatt sagt") print(" er, was eine zusaetzliche Montagestunde wert waere. Im Depot sagt") print(" dieselbe Zahl, was eine zusaetzliche Einheit Risikobudget wert") print(" waere - der PREIS DES RISIKOS. Das ist kein Sprachbild, sondern") print(" derselbe Dualwert derselben Nebenbedingung.") # --- Teil 2 ---------------------------------------------------------- print("\n" + "-" * 80) print("2. Absicherung gegen den schlechtesten Fall - EINE CVaR-Funktion\n") kurse, rendite, anlagen = kursszenarien(rng) anteile_depot = berichte_cvar( "Depot: die schlechtesten 5 % der Handelstage", kurse, rendite, anlagen, 5.5, "% je Tag ", "Erwartete Jahresrendite (%)") print() lieferung, zuverlaessig, lieferanten = lieferszenarien(rng) anteile_einkauf = berichte_cvar("Einkauf: die schlechtesten 5 % der Bestellungen", lieferung, zuverlaessig, lieferanten, 0.945, "Tage Verzug ", "Mittlere Zuverlaessigkeit") print("\n Auch hier: eine Funktion, zwei Aufrufe. Die Zielfunktion") print(" interessiert sich nicht dafuer, ob in der Matrix Prozentpunkte") print(" oder Tage stehen.") print() print(f" Interessant ist, WO die Streuungsgrenze von {MAX_ANTEIL:.0%} bindet:") print(f" Depot : {(anteile_depot > MAX_ANTEIL - 1e-4).sum()} von " f"{len(anteile_depot)} Posten am Anschlag") print(f" Einkauf : {(anteile_einkauf > MAX_ANTEIL - 1e-4).sum()} von " f"{len(anteile_einkauf)} Posten am Anschlag") print() print(" Im Einkauf zieht es die Loesung an den sicheren Lokalzulieferer,") print(" bis die Grenze sie stoppt. Im Depot nicht - dort verhindert die") print(" Mindestrendite, dass alles in Anleihen wandert. Zwei verschiedene") print(" Bremsen also, und sie stehen an verschiedenen Stellen des Modells:") print(" einmal in einer Nebenbedingung ueber die Anteile, einmal in einer") print(" ueber den Ertrag. Wer eine Struktur uebertraegt, uebertraegt eben") print(" nicht automatisch mit, WELCHE Bedingung am Ende bindet.") # --- Teil 3: Der ehrliche Teil --------------------------------------- print("\n" + "=" * 80) print(" WAS NICHT HINUEBERREICHT") print("=" * 80) print("Die Struktur traegt. Drei Unterschiede tragen NICHT mit, und wer sie") print("uebersieht, macht aus einer nuetzlichen Analogie einen Fehler:") print() print("1. WOHER DIE ZAHLEN KOMMEN.") print(" In der Werkstatt ist der Verbrauch je Tisch gemessen - drei") print(" Montagestunden sind drei Montagestunden. Im Depot ist die") print(" erwartete Rendite GESCHAETZT, und zwar mit einem Fehler, der") print(" groesser sein kann als die Unterschiede zwischen den Anlagen") print(" (Kapitel Markowitz, Renditeschaetzung_Falle.py). Dieselbe") print(" Optimierung ist im einen Fall Planung und im anderen") print(" Fehlerverstaerkung.") print() print("2. OB DIE VERGANGENHEIT ETWAS UEBER DIE ZUKUNFT SAGT.") print(" Lieferzeiten haben physikalische Ursachen: Entfernung, Zoll,") print(" Kapazitaet. Sie aendern sich langsam und nachvollziehbar.") print(" Kursrenditen entstehen aus dem Verhalten von Marktteilnehmern,") print(" die selbst auf Modelle reagieren - dort verschwindet ein") print(" erkanntes Muster oft genau deshalb, weil es erkannt wurde") print(" (Kapitel Handelsmaschine, Data_Snooping.py).") print() print("3. OB TEILBARKEIT ERLAUBT IST.") print(" 37,4 % eines Aktienfonds sind ein normaler Auftrag. 37,4 % eines") print(" Lieferanten sind es nicht - Vertraege, Mindestabnahmen und") print(" Ruestzeiten machen Einkaufsentscheidungen ganzzahlig. Genau") print(" deshalb ist Teil II voller MILP und Teil V fast frei davon.") print() print("Die Bruecke traegt also die MODELLE, nicht die Annahmen. Wer sie") print("benutzt, spart sich das Lernen der Methoden - nicht das Nachdenken") print("ueber die Daten.") print("=" * 80)