#!/usr/bin/env python3

# Entropie_Maximierte_Allokation.py
"""
Kapitel QP/NLP: Nichtlineare Optimierung (NLP) mit scipy SLSQP.
Modell: Mean-Variance-Portfolio mit Shannon-Entropie-Diversifikation.

Achtung, haeufiger Fehler: Sigma per A@A.T zu erzeugen und anschliessend die
Diagonale zu ueberschreiben zerstoert die positive Semidefinitheit - die
Matrix kann dadurch einen negativen Eigenwert bekommen und unmoegliche
Korrelationen (> 1) implizieren. Dieses Programm baut Sigma stattdessen als
D * C * D aus Volatilitaeten und einer echten Korrelationsmatrix -
konstruktionsbedingt immer PSD.

Zusaetzlich: Multistart, weil SLSQP nur lokale Optima findet.
"""

import numpy as np
import pandas as pd
from scipy.optimize import minimize

ASSETS = ["Tech-Aktien", "Rohstoffe", "US-Treasuries", "Krypto"]
MU = np.array([0.14, 0.07, 0.03, 0.22])          # erwartete Jahresrenditen
VOLA = np.array([0.20, 0.14, 0.07, 0.40])        # Volatilitaeten p.a.
KORRELATION = np.array([
    [1.00, 0.25, -0.10, 0.55],
    [0.25, 1.00, 0.05, 0.20],
    [-0.10, 0.05, 1.00, -0.15],
    [0.55, 0.20, -0.15, 1.00],
])

ALPHA = 1.0        # Gewicht des Ertrags
BETA = 0.015       # Gewicht der Entropie-Diversifikation
UNTERGRENZE = 0.001
OBERGRENZE = 0.80


def baue_kovarianz() -> np.ndarray:
    """Sigma = D * C * D. Immer PSD, wenn C eine gueltige Korrelationsmatrix ist."""
    eig_c = np.linalg.eigvalsh(KORRELATION)
    if eig_c.min() < -1e-10:
        raise ValueError(f"Korrelationsmatrix unmoeglich (Eigenwert {eig_c.min():.4f}).")
    D = np.diag(VOLA)
    Sigma = D @ KORRELATION @ D
    eig_s = np.linalg.eigvalsh(Sigma)
    assert eig_s.min() > 0, "Sigma nicht positiv definit!"
    return Sigma


SIGMA = baue_kovarianz()


def zielfunktion(w, alpha=ALPHA, beta=BETA):
    """Risiko - Ertrag - Entropiepraemie (wird minimiert)."""
    varianz = 0.5 * w @ SIGMA @ w
    ertrag = MU @ w
    entropie = -np.sum(w * np.log(w + 1e-12))
    return varianz - alpha * ertrag - beta * entropie


def gradient(w, alpha=ALPHA, beta=BETA):
    """Exakter analytischer Gradient - beschleunigt und stabilisiert SLSQP."""
    grad_varianz = SIGMA @ w                       # d/dw von 0.5 w'Sigma w
    grad_ertrag = -alpha * MU
    grad_entropie = beta * (np.log(w + 1e-12) + 1.0)
    return grad_varianz + grad_ertrag + grad_entropie


def optimiere(startpunkt, alpha=ALPHA, beta=BETA):
    """alpha und beta werden durchgereicht - so bleibt die Funktion seiteneffektfrei."""
    return minimize(
        zielfunktion, startpunkt, args=(alpha, beta), jac=gradient,
        method="SLSQP", bounds=[(UNTERGRENZE, OBERGRENZE)] * len(MU),
        constraints=({"type": "eq",
                      "fun": lambda w: np.sum(w) - 1.0,
                      "jac": lambda w: np.ones(len(MU))}),
        options={"ftol": 1e-12, "maxiter": 300})


if __name__ == "__main__":
    n = len(MU)
    eig = np.linalg.eigvalsh(SIGMA)

    print("=" * 74)
    print("   NICHTLINEARE ENTROPIE-OPTIMIERTE ASSET-ALLOKATION")
    print("=" * 74)
    print("Pruefung der Kovarianzmatrix:")
    print(f"  Eigenwerte:          {np.round(eig, 6)}")
    print(f"  Positiv definit:     {bool(eig.min() > 0)}")
    print(f"  Groesste Korrelation ausserhalb der Diagonale: "
          f"{np.abs(KORRELATION - np.eye(n)).max():.2f}  (muss <= 1 sein)")

    # --- Multistart: SLSQP findet nur lokale Optima -----------------------
    rng = np.random.default_rng(7)
    startpunkte = [np.ones(n) / n]                                  # Gleichgewichtung
    for _ in range(9):
        z = rng.random(n) + 0.05
        startpunkte.append(z / z.sum())

    ergebnisse = [optimiere(s) for s in startpunkte]
    erfolgreich = [r for r in ergebnisse if r.success]
    bestes = min(erfolgreich, key=lambda r: r.fun)
    zielwerte = np.array([r.fun for r in erfolgreich])

    print(f"\nMultistart mit {len(startpunkte)} Startpunkten:")
    print(f"  Erfolgreich konvergiert: {len(erfolgreich)}")
    print(f"  Spannweite der Zielwerte: {zielwerte.max() - zielwerte.min():.2e}")
    print("  -> " + ("alle Startpunkte fuehren zum selben Optimum (Indiz fuer Konvexitaet)"
                     if zielwerte.max() - zielwerte.min() < 1e-6
                     else "ACHTUNG: verschiedene lokale Optima gefunden!"))

    w_opt = bestes.x
    rendite = MU @ w_opt
    volatilitaet = np.sqrt(w_opt @ SIGMA @ w_opt)
    entropie = -np.sum(w_opt * np.log(w_opt))

    print(f"\nKonvergenz:                {bestes.message}")
    print(f"Iterationen:               {bestes.nit}")
    print(f"Erwartete Jahresrendite:   {rendite * 100:6.2f} %")
    print(f"Erwartete Volatilitaet:    {volatilitaet * 100:6.2f} %")
    print(f"Diversifikations-Entropie: {entropie:.4f} "
          f"(Maximum bei Gleichgewichtung: {np.log(n):.4f})")

    print("\n" + pd.DataFrame({
        "Asset": ASSETS,
        "Erw. Rendite": [f"{r*100:.1f} %" for r in MU],
        "Volatilitaet": [f"{v*100:.1f} %" for v in VOLA],
        "Gewicht": [f"{w*100:6.2f} %" for w in w_opt],
    }).to_string(index=False))

    # --- Vergleich: was passiert ohne Entropieterm? ----------------------
    ohne = optimiere(np.ones(n) / n, beta=0.0)
    print("\n--- Wirkung des Entropieterms ---")
    print(f"  Mit Entropie  (beta={BETA}): Gewichte {np.round(w_opt * 100, 1)}")
    print(f"  Ohne Entropie (beta=0):     Gewichte {np.round(ohne.x * 100, 1)}")
    print("  Der Entropieterm zieht Kapital aus der Spitzenposition heraus,")
    print("  ohne dass eine harte Obergrenze noetig waere.")
    print("=" * 74)