AI-Driven Risk Management for Nordic Power Futures and GoO Portfolios

import pandas as pd
import numpy as np

def create_temporal_features(df):
    """
    Create temporal features from datetime index.

    Args:
        df: DataFrame with DatetimeIndex

    Returns:
        DataFrame with additional temporal features
    """
    df = df.copy()

    # Basic calendar features
    df['hour'] = df.index.hour
    df['day_of_week'] = df.index.dayofweek
    df['month'] = df.index.month
    df['quarter'] = df.index.quarter

    # Cyclical encoding (important for avoiding discontinuity at boundaries)
    # Hour: 0-23 → sin/cos representation
    df['hour_sin'] = np.sin(2 * np.pi * df['hour'] / 24)
    df['hour_cos'] = np.cos(2 * np.pi * df['hour'] / 24)

    # Day of week: 0-6
    df['dow_sin'] = np.sin(2 * np.pi * df['day_of_week'] / 7)
    df['dow_cos'] = np.cos(2 * np.pi * df['day_of_week'] / 7)

    # Month: 1-12
    df['month_sin'] = np.sin(2 * np.pi * (df['month'] - 1) / 12)
    df['month_cos'] = np.cos(2 * np.pi * (df['month'] - 1) / 12)

    # Nordic-specific: Heating/cooling season
    # Northern latitudes have strong seasonal demand
    df['is_winter'] = df['month'].isin([11, 12, 1, 2, 3]).astype(int)
    df['is_summer'] = df['month'].isin([6, 7, 8]).astype(int)

    # Public holidays (reduces demand)
    # Simplified - in production use holiday package
    df['is_holiday'] = 0

    return df

def create_weather_features(df, temperature, wind_forecast, precipitation):
    """
    Create weather-derived features for power price modeling.

    Args:
        df: Base DataFrame with DatetimeIndex
        temperature: Series of temperature (°C)
        wind_forecast: Series of forecasted wind generation (MW)
        precipitation: Series of precipitation (mm/day)

    Returns:
        DataFrame with weather features
    """
    df = df.copy()

    # Temperature features
    df['temperature'] = temperature

    # Heating Degree Days (HDD) - base 18°C
    df['hdd'] = np.maximum(18 - temperature, 0)

    # Cooling Degree Days (CDD) - base 22°C
    df['cdd'] = np.maximum(temperature - 22, 0)

    # Temperature squared (capture non-linear effects)
    df['temp_squared'] = temperature ** 2

    # Temperature below freezing (triggers heating demand surge)
    df['temp_below_zero'] = (temperature < 0).astype(int)

    # Seasonal temperature deviation
    # Calculate rolling seasonal mean (365-day window)
    seasonal_mean = temperature.rolling(window=365, center=True).mean()
    df['temp_deviation'] = temperature - seasonal_mean

    # Wind features
    df['wind_forecast'] = wind_forecast

    # Wind capacity factor (0-1)
    installed_capacity = 10000  # MW - example value
    df['wind_capacity_factor'] = wind_forecast / installed_capacity

    # Wind forecast uncertainty (use ensemble spread if available)
    # Placeholder: use historical volatility as proxy
    df['wind_volatility'] = wind_forecast.rolling(window=168).std()  # 1 week

    # Precipitation features
    df['precipitation'] = precipitation

    # Cumulative precipitation (predicts hydro inflows)
    # 30-day rolling sum
    df['precip_30d'] = precipitation.rolling(window=30).sum()

    # 90-day rolling sum (seasonal reservoir filling)
    df['precip_90d'] = precipitation.rolling(window=90).sum()

    # Precipitation deviation from seasonal normal
    seasonal_precip = precipitation.rolling(window=365, center=True).mean()
    df['precip_deviation'] = precipitation - seasonal_precip

    return df

def create_fundamental_features(df, reservoir_level, gas_price, carbon_price, nuclear_capacity):
    """
    Create features from power market fundamentals.

    Args:
        df: Base DataFrame
        reservoir_level: TWh of stored hydro energy
        gas_price: EUR/MWh natural gas price
        carbon_price: EUR/ton CO2
        nuclear_capacity: MW available nuclear capacity

    Returns:
        DataFrame with fundamental features
    """
    df = df.copy()

    # Reservoir levels
    max_capacity = 120  # TWh for Nordic region
    df['reservoir_twh'] = reservoir_level
    df['reservoir_filling'] = reservoir_level / max_capacity

    # Reservoir filling relative to historical average for this time of year
    seasonal_avg = reservoir_level.groupby([df.index.month, df.index.day]).transform('mean')
    df['reservoir_deviation'] = reservoir_level - seasonal_avg

    # Low reservoir indicator (critical level)
    df['reservoir_critical'] = (df['reservoir_filling'] < 0.3).astype(int)

    # Rate of reservoir change (weekly change)
    df['reservoir_change_weekly'] = reservoir_level - reservoir_level.shift(168)  # 1 week

    # Fuel prices
    df['gas_price'] = gas_price
    df['carbon_price'] = carbon_price

    # Thermal generation cost (gas + carbon)
    # Simplified: Gas plant efficiency ~50%, carbon intensity ~0.4 ton/MWh
    thermal_cost = gas_price / 0.5 + carbon_price * 0.4
    df['thermal_generation_cost'] = thermal_cost

    # Nuclear availability
    df['nuclear_capacity'] = nuclear_capacity
    df['nuclear_capacity_pct'] = nuclear_capacity / 10000  # % of installed capacity

    return df

def create_transmission_features(df, no_se_flow, no_se_capacity, no2_price, no5_price):
    """
    Create transmission/interconnection features.

    Args:
        df: Base DataFrame
        no_se_flow: MW flow from Norway to Sweden (positive = export)
        no_se_capacity: MW maximum capacity
        no2_price: Price in NO2 (Southern Norway)
        no5_price: Price in NO5 (Western Norway)

    Returns:
        DataFrame with transmission features
    """
    df = df.copy()

    # Interconnector utilization
    df['no_se_flow'] = no_se_flow
    df['no_se_utilization'] = np.abs(no_se_flow) / no_se_capacity

    # Congestion indicator (>90% utilization)
    df['no_se_congested'] = (df['no_se_utilization'] > 0.9).astype(int)

    # Price spreads (basis risk)
    df['no2_price'] = no2_price
    df['no5_price'] = no5_price
    df['no2_no5_spread'] = no2_price - no5_price

    # Spread volatility (indicates transmission risk)
    df['spread_volatility'] = df['no2_no5_spread'].rolling(window=168).std()

    return df

import pandas as pd
import numpy as np
from lightgbm import LGBMRegressor
from sklearn.model_selection import TimeSeriesSplit
from sklearn.metrics import mean_absolute_error, mean_squared_error

def train_price_forecast_model(df, target_col='price', n_lags=24):
    """
    Train gradient boosting model for power price forecasting.

    Args:
        df: DataFrame with features and target
        target_col: Name of target column (price)
        n_lags: Number of lagged price features to include

    Returns:
        Trained model and feature names
    """
    # Create lagged price features
    for lag in range(1, n_lags + 1):
        df[f'price_lag_{lag}'] = df[target_col].shift(lag)

    # Create rolling statistics
    df['price_rolling_mean_24h'] = df[target_col].shift(1).rolling(window=24).mean()
    df['price_rolling_std_24h'] = df[target_col].shift(1).rolling(window=24).std()
    df['price_rolling_mean_168h'] = df[target_col].shift(1).rolling(window=168).mean()

    # Drop NaN rows from lagging
    df = df.dropna()

    # Split features and target
    feature_cols = [col for col in df.columns if col != target_col and col != 'price']
    X = df[feature_cols]
    y = df[target_col]

    # Time series split for validation
    tscv = TimeSeriesSplit(n_splits=5)

    # Train model
    model = LGBMRegressor(
        n_estimators=500,
        learning_rate=0.05,
        max_depth=8,
        num_leaves=64,
        min_child_samples=20,
        subsample=0.8,
        colsample_bytree=0.8,
        reg_alpha=0.1,  # L1 regularization
        reg_lambda=1.0,  # L2 regularization
        random_state=42,
        n_jobs=-1
    )

    # Cross-validation to assess performance
    cv_scores = []
    for train_idx, val_idx in tscv.split(X):
        X_train, X_val = X.iloc[train_idx], X.iloc[val_idx]
        y_train, y_val = y.iloc[train_idx], y.iloc[val_idx]

        model.fit(
            X_train, y_train,
            eval_set=[(X_val, y_val)],
            eval_metric='mae',
            early_stopping_rounds=50,
            verbose=False
        )

        y_pred = model.predict(X_val)
        mae = mean_absolute_error(y_val, y_pred)
        cv_scores.append(mae)

    print(f"Cross-validation MAE: {np.mean(cv_scores):.2f} ± {np.std(cv_scores):.2f} EUR/MWh")

    # Train final model on all data
    model.fit(X, y, verbose=False)

    # Feature importance
    feature_importance = pd.DataFrame({
        'feature': feature_cols,
        'importance': model.feature_importances_
    }).sort_values('importance', ascending=False)

    print("\nTop 10 most important features:")
    print(feature_importance.head(10))

    return model, feature_cols

# Example usage:
# Assuming df contains all engineered features
# model, features = train_price_forecast_model(df, target_col='system_price')

import numpy as np
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers
from sklearn.preprocessing import StandardScaler

def create_lstm_dataset(data, target, lookback=168, horizon=24):
    """
    Create sequences for LSTM training.

    Args:
        data: Feature array (n_samples, n_features)
        target: Target array (n_samples,)
        lookback: Number of past timesteps to use
        horizon: Number of future timesteps to predict

    Returns:
        X: (n_samples, lookback, n_features)
        y: (n_samples, horizon)
    """
    X, y = [], []

    for i in range(lookback, len(data) - horizon + 1):
        X.append(data[i-lookback:i])
        y.append(target[i:i+horizon])

    return np.array(X), np.array(y)

def build_lstm_model(input_shape, horizon=24):
    """
    Build LSTM model for multi-step price forecasting.

    Args:
        input_shape: (lookback, n_features)
        horizon: Number of steps to forecast

    Returns:
        Compiled Keras model
    """
    model = keras.Sequential([
        # First LSTM layer with return sequences
        layers.LSTM(128, return_sequences=True, input_shape=input_shape),
        layers.Dropout(0.2),

        # Second LSTM layer
        layers.LSTM(64, return_sequences=False),
        layers.Dropout(0.2),

        # Dense layers
        layers.Dense(64, activation='relu'),
        layers.Dropout(0.1),

        # Output layer (forecasting 'horizon' steps)
        layers.Dense(horizon)
    ])

    model.compile(
        optimizer=keras.optimizers.Adam(learning_rate=0.001),
        loss='mse',
        metrics=['mae']
    )

    return model

def train_lstm_forecaster(df, feature_cols, target_col='price',
                          lookback=168, horizon=24):
    """
    Train LSTM model for power price forecasting.

    Args:
        df: DataFrame with features and target
        feature_cols: List of feature column names
        target_col: Target column name
        lookback: Hours of history to use (168 = 1 week)
        horizon: Hours ahead to forecast (24 = 1 day)

    Returns:
        Trained model and scalers
    """
    # Prepare data
    feature_data = df[feature_cols].values
    target_data = df[target_col].values

    # Scale features (important for neural networks)
    feature_scaler = StandardScaler()
    feature_data_scaled = feature_scaler.fit_transform(feature_data)

    target_scaler = StandardScaler()
    target_data_scaled = target_scaler.fit_transform(target_data.reshape(-1, 1)).flatten()

    # Create sequences
    X, y = create_lstm_dataset(feature_data_scaled, target_data_scaled,
                                lookback=lookback, horizon=horizon)

    # Split train/validation (80/20, maintaining temporal order)
    split_idx = int(0.8 * len(X))
    X_train, X_val = X[:split_idx], X[split_idx:]
    y_train, y_val = y[:split_idx], y[split_idx:]

    print(f"Training set: {X_train.shape}, Validation set: {X_val.shape}")

    # Build and train model
    model = build_lstm_model(input_shape=(lookback, len(feature_cols)), horizon=horizon)

    early_stopping = keras.callbacks.EarlyStopping(
        monitor='val_loss',
        patience=10,
        restore_best_weights=True
    )

    reduce_lr = keras.callbacks.ReduceLROnPlateau(
        monitor='val_loss',
        factor=0.5,
        patience=5,
        min_lr=1e-6
    )

    history = model.fit(
        X_train, y_train,
        validation_data=(X_val, y_val),
        epochs=100,
        batch_size=64,
        callbacks=[early_stopping, reduce_lr],
        verbose=1
    )

    # Evaluate
    val_pred_scaled = model.predict(X_val)

    # Inverse transform for multi-horizon predictions
    # Reshape to (n_samples * horizon, 1), inverse transform, reshape back
    n_samples = val_pred_scaled.shape[0]
    val_pred = target_scaler.inverse_transform(
        val_pred_scaled.reshape(-1, 1)
    ).reshape(n_samples, horizon)

    val_true = target_scaler.inverse_transform(
        y_val.reshape(-1, 1)
    ).reshape(n_samples, horizon)

    mae = np.mean(np.abs(val_pred - val_true))
    print(f"\nValidation MAE (multi-step): {mae:.2f} EUR/MWh")

    return model, feature_scaler, target_scaler

# Example usage:
# model, feat_scaler, targ_scaler = train_lstm_forecaster(
#     df, feature_cols=['temperature', 'wind_forecast', 'reservoir_twh', ...],
#     target_col='system_price'
# )

import lightgbm as lgb

def train_quantile_models(df, feature_cols, target_col='price',
                           quantiles=[0.01, 0.05, 0.25, 0.50, 0.75, 0.95, 0.99]):
    """
    Train multiple quantile regression models to estimate price distribution.

    Args:
        df: DataFrame with features
        feature_cols: Feature column names
        target_col: Target column name
        quantiles: List of quantiles to predict

    Returns:
        Dictionary of models, one per quantile
    """
    X = df[feature_cols]
    y = df[target_col]

    # Time series split
    split_idx = int(0.8 * len(X))
    X_train, X_val = X.iloc[:split_idx], X.iloc[split_idx:]
    y_train, y_val = y.iloc[:split_idx], y.iloc[split_idx:]

    models = {}

    for q in quantiles:
        print(f"\nTraining quantile {q:.2f} model...")

        model = lgb.LGBMRegressor(
            objective='quantile',
            alpha=q,  # Quantile to predict
            n_estimators=300,
            learning_rate=0.05,
            max_depth=8,
            num_leaves=64,
            random_state=42
        )

        model.fit(
            X_train, y_train,
            eval_set=[(X_val, y_val)],
            eval_metric='quantile',
            early_stopping_rounds=30,
            verbose=False
        )

        models[q] = model

        # Evaluate coverage (% of actual values below predicted quantile)
        val_pred = model.predict(X_val)
        coverage = (y_val <= val_pred).mean()
        print(f"  Quantile {q:.2f}: Coverage = {coverage:.3f} (target = {q:.2f})")

    return models

def predict_price_distribution(models, X_new):
    """
    Predict full price distribution for new data.

    Args:
        models: Dictionary of quantile models
        X_new: New feature data

    Returns:
        DataFrame with quantile predictions
    """
    predictions = {}

    for q, model in models.items():
        predictions[f'q{int(q*100)}'] = model.predict(X_new)

    return pd.DataFrame(predictions)

# Example usage:
# quantile_models = train_quantile_models(df, feature_cols)
# future_distribution = predict_price_distribution(quantile_models, X_test)

def calculate_portfolio_var(quantile_predictions, positions, confidence=0.95):
    """
    Calculate Value-at-Risk for power portfolio using quantile forecasts.

    Args:
        quantile_predictions: DataFrame with columns ['q01', 'q05', ..., 'q99']
                              rows = hours/products
        positions: Array of position sizes (MW) matching rows of quantile_predictions
        confidence: VaR confidence level (0.95 = 95%)

    Returns:
        VaR value (positive = potential loss)
    """
    # Select appropriate quantile column
    # For 95% VaR, use 5th percentile (0.05 quantile)
    alpha = 1 - confidence
    quantile_col = f'q{int(alpha * 100)}'

    if quantile_col not in quantile_predictions.columns:
        raise ValueError(f"Quantile {quantile_col} not available")

    # Calculate portfolio value at quantile
    prices_at_quantile = quantile_predictions[quantile_col].values
    portfolio_value_at_quantile = np.sum(positions * prices_at_quantile)

    # Calculate expected value
    expected_prices = quantile_predictions['q50'].values  # Median
    expected_portfolio_value = np.sum(positions * expected_prices)

    # VaR is the loss relative to expected value
    var = expected_portfolio_value - portfolio_value_at_quantile

    return var

# Example:
# positions = np.array([100, 150, -50, ...])  # MW positions (negative = short)
# var_95 = calculate_portfolio_var(quantile_predictions, positions, confidence=0.95)
# print(f"95% VaR: €{var_95:,.0f}")

def calculate_portfolio_cvar(quantile_predictions, positions, confidence=0.95):
    """
    Calculate Conditional Value-at-Risk (CVaR) using quantile forecasts.

    Args:
        quantile_predictions: DataFrame with quantile predictions
        positions: Portfolio positions (MW)
        confidence: CVaR confidence level

    Returns:
        CVaR value (expected loss in tail)
    """
    alpha = 1 - confidence

    # Get all quantiles in the tail (q <= alpha)
    tail_quantiles = [col for col in quantile_predictions.columns
                      if col.startswith('q') and int(col[1:]) <= int(alpha * 100)]

    if not tail_quantiles:
        raise ValueError(f"No quantiles available for tail probability {alpha}")

    # Calculate portfolio values at each tail quantile
    tail_portfolio_values = []
    for q_col in tail_quantiles:
        prices = quantile_predictions[q_col].values
        portfolio_value = np.sum(positions * prices)
        tail_portfolio_values.append(portfolio_value)

    # CVaR is mean portfolio value in tail
    cvar_portfolio_value = np.mean(tail_portfolio_values)

    # Express as loss relative to expected value
    expected_value = np.sum(positions * quantile_predictions['q50'].values)
    cvar = expected_value - cvar_portfolio_value

    return cvar

# Example:
# cvar_95 = calculate_portfolio_cvar(quantile_predictions, positions, confidence=0.95)
# print(f"95% CVaR: €{cvar_95:,.0f}")
# print(f"CVaR exceeds VaR by: €{(cvar_95 - var_95):,.0f}")

def monte_carlo_var_cvar(forecast_model, feature_data, positions,
                          n_scenarios=10000, confidence=0.95):
    """
    Calculate VaR and CVaR using Monte Carlo simulation.

    Args:
        forecast_model: ML model that predicts price mean and volatility
        feature_data: Features for forecasting period
        positions: Portfolio positions (MW)
        n_scenarios: Number of scenarios to simulate
        confidence: Risk confidence level

    Returns:
        var, cvar values
    """
    # Get point forecast and uncertainty from model
    # Assuming model predicts both mean and std
    price_mean = forecast_model['mean'].predict(feature_data)
    price_std = forecast_model['std'].predict(feature_data)

    n_hours = len(price_mean)

    # Build covariance matrix with serial correlation
    # Assume AR(1) structure: corr(t, t+k) = rho^k
    # Calibrate rho from historical residual ACF; 0.7 is illustrative only
    rho = 0.7  # <--- replace with estimated autocorrelation from residuals

    # Covariance matrix with exponential decay
    time_diff = np.abs(np.arange(n_hours)[:, None] - np.arange(n_hours)[None, :])
    correlation_matrix = rho ** time_diff

    # Convert correlation to covariance
    # price_std should be marginal hourly volatility (not a wide CI band)
    std_matrix = np.outer(price_std, price_std)
    covariance_matrix = correlation_matrix * std_matrix

    # Generate correlated scenarios using multivariate normal
    # NOTE: Using normal (not log-normal) to allow negative prices
    # Nordic markets do experience negative prices during high wind/low demand
    price_scenarios = np.random.multivariate_normal(
        mean=price_mean,
        cov=covariance_matrix,
        size=n_scenarios
    )

    # Calculate portfolio value for each scenario
    portfolio_values = price_scenarios @ positions  # Matrix multiplication

    # Calculate risk metrics
    portfolio_pnl = portfolio_values - np.mean(portfolio_values)

    # VaR: loss at (1-confidence) percentile
    var = -np.percentile(portfolio_pnl, (1 - confidence) * 100)

    # CVaR: mean loss beyond VaR
    tail_losses = -portfolio_pnl[portfolio_pnl < -var]
    cvar = np.mean(tail_losses) if len(tail_losses) > 0 else var

    return var, cvar, price_scenarios

# Example:
# var, cvar, scenarios = monte_carlo_var_cvar(
#     forecast_model, X_tomorrow, positions, n_scenarios=10000
# )
# print(f"Monte Carlo VaR: €{var:,.0f}")
# print(f"Monte Carlo CVaR: €{cvar:,.0f}")

import numpy as np
from scipy.optimize import minimize

class HydroPortfolioManager:
    """
    AI-driven portfolio manager for hydro power producer.
    """

    def __init__(self, capacity_mw, reservoir_capacity_twh,
                 inflow_model, price_model):
        self.capacity_mw = capacity_mw
        self.reservoir_capacity_twh = reservoir_capacity_twh
        self.inflow_model = inflow_model
        self.price_model = price_model

    def forecast_generation(self, current_month, n_months=12):
        """
        Forecast generation for next n months based on inflow model.

        Returns:
            Array of forecasted monthly generation (MWh)
        """
        # Get inflow forecasts from ML model
        inflow_forecasts = self.inflow_model.predict(current_month, n_months)

        # Convert inflows to generation (accounting for reservoir dynamics)
        generation = np.zeros(n_months)
        reservoir = self.reservoir_capacity_twh * 1000  # MWh

        for month in range(n_months):
            # Inflow adds to reservoir
            reservoir += inflow_forecasts[month]

            # Generate up to capacity or available water
            hours_in_month = 730  # Average
            max_generation = self.capacity_mw * hours_in_month
            available_generation = min(max_generation, reservoir)

            generation[month] = available_generation
            reservoir -= generation[month]

            # Reservoir can't exceed capacity
            reservoir = min(reservoir, self.reservoir_capacity_twh * 1000)

        return generation

    def forecast_price_distribution(self, generation_forecast, current_features):
        """
        Forecast price distribution given generation forecast.

        Returns:
            Dictionary with quantile forecasts
        """
        # Use ML price model to forecast distribution
        # Input features include forecasted generation (supply proxy)
        features = current_features.copy()
        features['forecasted_generation'] = generation_forecast.mean()

        price_quantiles = {}
        for q in [0.05, 0.25, 0.50, 0.75, 0.95]:
            price_quantiles[q] = self.price_model[q].predict([features])[0]

        return price_quantiles

    def optimize_hedge_ratio(self, generation_forecast, price_distribution,
                             forward_price, risk_aversion=2.0):
        """
        Optimize hedge ratio to maximize risk-adjusted return.

        Args:
            generation_forecast: Array of forecasted generation (MWh)
            price_distribution: Dict with price quantiles
            forward_price: Current forward price for hedging (EUR/MWh)
            risk_aversion: Coefficient for CVaR penalty (higher = more conservative)

        Returns:
            Optimal hedge ratio (0-1)
        """
        def objective(hedge_ratio):
            # Portfolio value = spot revenue + hedge P&L

            # Spot revenue under different price scenarios
            spot_scenarios = np.array([
                price_distribution[0.05],
                price_distribution[0.25],
                price_distribution[0.50],
                price_distribution[0.75],
                price_distribution[0.95]
            ])

            generation = generation_forecast.sum()
            spot_revenue = spot_scenarios * generation

            # Hedge P&L: gain if spot < forward, loss if spot > forward
            hedge_volume = hedge_ratio * generation
            hedge_pnl = (forward_price - spot_scenarios) * hedge_volume

            # Total portfolio value
            portfolio_value = spot_revenue + hedge_pnl

            # Risk-adjusted objective: Expected value - risk_aversion * CVaR
            expected_value = np.mean(portfolio_value)

            # CVaR: mean of worst 20% scenarios
            # With 5 scenarios, worst 20% = 1 scenario (the worst one).
            # For production, increase scenarios (e.g., 100+) to make CVaR smoother.
            worst_scenarios = np.partition(portfolio_value, 0)[:1]
            cvar = -np.mean(worst_scenarios)  # Negative because we want loss

            return -(expected_value - risk_aversion * cvar)  # Minimize negative

        # Optimize hedge ratio between 0 and 1
        result = minimize(
            objective,
            x0=[0.7],  # Initial guess
            bounds=[(0, 1)],
            method='L-BFGS-B'
        )

        return result.x[0]

    def monthly_rebalance(self, current_month, current_features, forward_price):
        """
        Monthly portfolio rebalancing decision.

        Returns:
            Recommended hedge ratio and risk metrics
        """
        # Forecast generation
        generation_forecast = self.forecast_generation(current_month)

        # Forecast prices
        price_dist = self.forecast_price_distribution(
            generation_forecast, current_features
        )

        # Optimize hedge
        optimal_hedge = self.optimize_hedge_ratio(
            generation_forecast, price_dist, forward_price
        )

        print(f"\n=== Monthly Portfolio Review ===")
        print(f"Forecasted generation: {generation_forecast.sum():,.0f} MWh")
        print(f"Price forecast (median): €{price_dist[0.50]:.2f}/MWh")
        print(f"Price forecast (5th-95th percentile): €{price_dist[0.05]:.2f} - €{price_dist[0.95]:.2f}/MWh")
        print(f"Current forward price: €{forward_price:.2f}/MWh")
        print(f"\nRecommended hedge ratio: {optimal_hedge:.1%}")

        return {
            'hedge_ratio': optimal_hedge,
            'generation_forecast': generation_forecast,
            'price_distribution': price_dist
        }

# Example usage:
# manager = HydroPortfolioManager(
#     capacity_mw=500,
#     reservoir_capacity_twh=3.0,
#     inflow_model=lstm_inflow_model,
#     price_model=quantile_price_models
# )
#
# decision = manager.monthly_rebalance(
#     current_month=pd.Timestamp('2024-01-01'),
#     current_features=current_feature_vector,
#     forward_price=45.0
# )

from sklearn.mixture import GaussianMixture
import matplotlib.pyplot as plt

class GoORegimeSwitchingModel:
    """
    Regime-switching model for GoO certificate pricing and risk.
    """

    def __init__(self, n_regimes=3):
        self.n_regimes = n_regimes
        self.gmm = None

    def fit(self, goo_prices, power_prices, renewable_capacity, corporate_demand):
        """
        Fit Gaussian Mixture Model to identify GoO market regimes.

        NOTE: This uses unsupervised clustering (GaussianMixture) rather than a true
        Hidden Markov Model. It identifies distinct market states but does NOT model
        temporal persistence or transition probabilities between regimes. For a true
        regime-switching model with temporal dynamics, use hmmlearn.GaussianHMM instead.
        The class name is kept for continuity; swap the estimator if you need an HMM.

        Args:
            goo_prices: Time series of GoO prices (EUR/MWh)
            power_prices: Time series of power prices
            renewable_capacity: Installed renewable capacity (MW)
            corporate_demand: Corporate GoO demand proxy (# of announcements, etc.)

        Returns:
            Fitted model with identified regimes
        """
        # Features for regime identification
        features = np.column_stack([
            goo_prices,
            power_prices,
            goo_prices / (power_prices + 1),  # GoO as % of power price
            renewable_capacity,
            corporate_demand,
            np.gradient(goo_prices),  # GoO price momentum
        ])

        # Fit Gaussian Mixture Model to identify regimes via clustering
        # This groups similar market conditions together but doesn't enforce
        # temporal smoothness (i.e., regimes can switch rapidly day-to-day)
        self.gmm = GaussianMixture(
            n_components=self.n_regimes,
            covariance_type='full',
            random_state=42
        )

        self.gmm.fit(features)

        # Predict regime for each time period
        regimes = self.gmm.predict(features)

        # Characterize each regime
        regime_stats = {}
        for r in range(self.n_regimes):
            regime_mask = (regimes == r)
            regime_stats[r] = {
                'mean_goo_price': goo_prices[regime_mask].mean(),
                'mean_power_price': power_prices[regime_mask].mean(),
                'goo_volatility': goo_prices[regime_mask].std(),
                'power_goo_correlation': np.corrcoef(
                    power_prices[regime_mask],
                    goo_prices[regime_mask]
                )[0, 1],
                'frequency': regime_mask.sum() / len(regimes)
            }

        self.regime_stats = regime_stats

        print("=== Identified GoO Market Regimes ===")
        for r, stats in regime_stats.items():
            print(f"\nRegime {r} ({stats['frequency']:.1%} of time):")
            print(f"  Mean GoO price: €{stats['mean_goo_price']:.2f}/MWh")
            print(f"  Mean power price: €{stats['mean_power_price']:.2f}/MWh")
            print(f"  GoO volatility: €{stats['goo_volatility']:.2f}/MWh")
            print(f"  Power-GoO correlation: {stats['power_goo_correlation']:.2f}")

        return self

    def predict_regime(self, current_features):
        """
        Predict current market regime.

        Args:
            current_features: Current values of regime-identifying features

        Returns:
            Most likely regime (0 to n_regimes-1)
        """
        return self.gmm.predict([current_features])[0]

    def regime_specific_var(self, portfolio_positions, current_regime):
        """
        Calculate VaR using regime-specific correlations.

        Args:
            portfolio_positions: Dict with 'power_mwh' and 'goo_certificates'
            current_regime: Current market regime

        Returns:
            VaR estimate
        """
        regime_stats = self.regime_stats[current_regime]

        # Use regime-specific volatilities and correlation
        power_vol = 15.0  # Simplified - should come from model
        goo_vol = regime_stats['goo_volatility']
        correlation = regime_stats['power_goo_correlation']

        # Portfolio variance
        power_var = (portfolio_positions['power_mwh'] * power_vol) ** 2
        goo_var = (portfolio_positions['goo_certificates'] * goo_vol) ** 2
        covar = 2 * portfolio_positions['power_mwh'] * portfolio_positions['goo_certificates'] * power_vol * goo_vol * correlation

        portfolio_std = np.sqrt(power_var + goo_var + covar)

        # VaR at 95% confidence (1.645 standard deviations)
        var_95 = 1.645 * portfolio_std

        return var_95

# Example usage:
# model = GoORegimeSwitchingModel(n_regimes=3)
# model.fit(goo_price_history, power_price_history, renewable_capacity_history, corporate_demand_history)
#
# current_regime = model.predict_regime(current_feature_vector)
# var = model.regime_specific_var(
#     {'power_mwh': 100000, 'goo_certificates': 100000},
#     current_regime
# )

def optimize_power_goo_portfolio(power_forecast, goo_forecast,
                                  power_forward, goo_forward,
                                  risk_aversion=1.5):
    """
    Optimize combined power and GoO hedging strategy.

    Args:
        power_forecast: ML forecast of power price distribution
        goo_forecast: ML forecast of GoO price distribution
        power_forward: Current power forward price
        goo_forward: Current GoO forward price
        risk_aversion: Risk aversion parameter

    Returns:
        Optimal hedge ratios for power and GoO
    """
    def objective(hedges):
        power_hedge, goo_hedge = hedges

        # Generate correlated scenarios
        n_scenarios = 1000
        # Sample from forecasted distributions (simplified - use copula in practice)
        power_scenarios = np.random.normal(
            power_forecast['mean'],
            power_forecast['std'],
            n_scenarios
        )
        goo_scenarios = np.random.normal(
            goo_forecast['mean'],
            goo_forecast['std'],
            n_scenarios
        )

        # Apply correlation
        correlation = goo_forecast.get('power_goo_corr', 0.3)
        goo_scenarios = (correlation * power_scenarios +
                         np.sqrt(1 - correlation**2) * goo_scenarios)

        # Portfolio value under each scenario
        generation = 10000  # MWh
        certificates = 10000  # Equal to generation

        spot_revenue = power_scenarios * generation + goo_scenarios * certificates
        hedge_pnl = ((power_forward - power_scenarios) * power_hedge * generation +
                     (goo_forward - goo_scenarios) * goo_hedge * certificates)

        total_value = spot_revenue + hedge_pnl

        # Risk-adjusted return
        expected = np.mean(total_value)
        cvar = -np.mean(np.partition(total_value, int(0.05 * n_scenarios))[:int(0.05 * n_scenarios)])

        return -(expected - risk_aversion * cvar)

    # Optimize
    from scipy.optimize import minimize
    result = minimize(
        objective,
        x0=[0.5, 0.5],
        bounds=[(0, 1), (0, 1)],
        method='L-BFGS-B'
    )

    return {'power_hedge': result.x[0], 'goo_hedge': result.x[1]}

import shap

def explain_risk_prediction(model, features, feature_names):
    """
    Generate SHAP explanations for risk model predictions.

    Args:
        model: Trained ML model
        features: Feature values for scenario to explain
        feature_names: Names of features

    Returns:
        SHAP values showing feature contributions
    """
    # Create SHAP explainer
    explainer = shap.TreeExplainer(model)  # For tree-based models

    # Calculate SHAP values
    shap_values = explainer.shap_values(features)

    # Visualize
    shap.summary_plot(shap_values, features, feature_names=feature_names)

    # Get top contributors
    feature_importance = np.abs(shap_values).mean(axis=0)
    top_features = pd.DataFrame({
        'feature': feature_names,
        'importance': feature_importance
    }).sort_values('importance', ascending=False)

    print("Top contributors to risk prediction:")
    print(top_features.head(10))

    return shap_values