Portefeuille Markowitz Covariance Dynamique

# Excerpt: official allocation rule. Trading costs are applied separately.
import math
from typing import Any
import numpy as np
from scipy.optimize import minimize
from sklearn.covariance import LedoitWolf
TICKERS = ("SPY", "EFA", "IEF", "VNQ", "GLD")
ANNUAL_RISK_FREE = 0.02  # optimization parameter only
TRANSACTION_COST_RATE = 0.0009

class ReplicationError(ValueError):
    pass

def optimize_weights(window_returns: np.ndarray, tickers: tuple[str, ...] = TICKERS) -> tuple[np.ndarray, dict[str, Any]]:
    """Optimize annualized plug-in Sharpe using deterministic feasible SLSQP starts."""
    matrix = np.asarray(window_returns, dtype=float)
    if matrix.shape != (252, len(tickers)) or not np.isfinite(matrix).all():
        raise ReplicationError("Optimization requires exactly 252 finite daily return rows")
    try:
        estimator = LedoitWolf().fit(matrix)
        daily_cov = np.asarray(estimator.covariance_, dtype=float)
        annual_cov = daily_cov * 252.0
        annual_mean = matrix.mean(axis=0) * 252.0
    except Exception as exc:
        raise ReplicationError(f"Ledoit-Wolf covariance estimation failed: {exc}") from exc
    if not np.isfinite(annual_mean).all() or not np.isfinite(annual_cov).all():
        raise ReplicationError("Non-finite expected return/covariance estimate")
    if not np.allclose(annual_cov, annual_cov.T, atol=1e-12):
        raise ReplicationError("Estimated covariance is not symmetric")
    eigenvalues = np.linalg.eigvalsh(annual_cov)
    if not np.isfinite(eigenvalues).all() or eigenvalues.min() <= 0:
        raise ReplicationError("Estimated covariance is not positive definite")

    def objective(weights: np.ndarray) -> float:
        variance = float(weights @ annual_cov @ weights)
        if not math.isfinite(variance) or variance <= 0:
            return 1e100
        value = (float(weights @ annual_mean) - ANNUAL_RISK_FREE) / math.sqrt(variance)
        return -value if math.isfinite(value) else 1e100

    n = len(tickers)
    starts = [np.full(n, 1.0 / n)]
    # Deterministic bounded vertices and two-asset mixes provide repeatable
    # multi-start coverage of the capped simplex without stochastic seeds.
    for high in range(n):
        candidate = np.full(n, 0.125)
        candidate[high] = 0.5
        starts.append(candidate)
    for left in range(n):
        for right in range(left + 1, n):
            candidate = np.zeros(n)
            candidate[left] = candidate[right] = 0.5
            starts.append(candidate)
    constraints = ({"type": "eq", "fun": lambda w: float(np.sum(w) - 1.0)},)
    feasible: list[tuple[float, np.ndarray, int]] = []
    solver_errors = []
    for start in starts:
        try:
            result = minimize(
                objective, start, method="SLSQP", bounds=[(0.0, 0.5)] * n,
                constraints=constraints,
                options={"ftol": 1e-12, "maxiter": 1000, "disp": False},
            )
        except Exception as exc:
            solver_errors.append(str(exc))
            continue
        weights = np.asarray(result.x, dtype=float)
        is_feasible = (
            weights.shape == (n,) and np.isfinite(weights).all()
            and abs(float(weights.sum()) - 1.0) <= 1e-7
            and (weights >= -1e-7).all() and (weights <= 0.5 + 1e-7).all()
        )
        if result.success and is_feasible and math.isfinite(objective(weights)):
            weights = np.clip(weights, 0.0, 0.5)
            weights /= weights.sum()
            feasible.append((objective(weights), weights, int(result.nit)))
        elif not result.success:
            solver_errors.append(str(result.message))
    if not feasible:
        raise ReplicationError("All deterministic SLSQP starts failed: " + "; ".join(solver_errors[:5]))
    feasible.sort(key=lambda item: (item[0], tuple(float(x) for x in item[1])))
    score, weights, iterations = feasible[0]
    # Numerical solver tolerances can leave a bound off by a few ulps. Project
    # onto the capped simplex without a normalize-after-clipping cap violation.
    low_shift, high_shift = -1.0, 1.0
    for _ in range(100):
        shift = (low_shift + high_shift) / 2
        projected = np.clip(weights + shift, 0.0, 0.5)
        if projected.sum() < 1.0:
            low_shift = shift
        else:
            high_shift = shift
    weights = np.clip(weights + (low_shift + high_shift) / 2, 0.0, 0.5)
    if abs(float(weights.sum()) - 1.0) > 1e-10:
        raise ReplicationError("Could not project SLSQP result onto the capped simplex")
    return weights, {
        "daily_return_observations": 252,
        "annual_expected_returns_arithmetic_mean": dict(zip(tickers, annual_mean.tolist())),
        "annual_covariance_estimator": "sklearn.covariance.LedoitWolf.fit(Tx5 daily simple returns) * 252",
        "optimized_annualized_sharpe": -score,
        "solver": "scipy.optimize.minimize(method=SLSQP)",
        "solver_multistarts": len(starts),
        "selected_solver_iterations": iterations,
    }