Adding CVaR Optimization using GARCH-LSTM-EVT

[Milan933-coder]

GARCH-EVT-LSTM Hybrid Simulator with CVaR Optimization

📌 Description

While Qlib offers a robust Model Zoo of deep learning architectures (LSTMs, Transformers, LightGBM), it currently lacks hybrid models that explicitly combine deep learning with rigorous stochastic volatility and extreme value statistical frameworks.

This feature proposes the integration of a GARCH-EVT-LSTM Hybrid Simulator coupled with a Conditional Value at Risk (CVaR) Portfolio Optimizer.

Standard deep learning models often fail to capture the fat tails and volatility clustering inherent in financial markets. By integrating Extreme Value Theory (EVT) and GARCH, this model stabilizes the LSTM's learning process and provides robust quantile forecasts, making it resilient to black-swan events (e.g., COVID-19 crash).

---

🚀 Proposed Core Components

1. EVT (Extreme Value Theory) Smoothing Module

• Uses Peak-Over-Threshold (POT)
• Fits Generalized Pareto Distribution (GPD) on extreme residuals
• Prevents overfitting to anomalies while preserving tail risk

2. LSTM-GARCH Joint Simulator

• GARCH(1,1) → models conditional heteroskedasticity (volatility clustering)
• LSTM → predicts conditional mean of returns
• Cholesky / Copula → models cross-asset dependencies

3. CVaR Portfolio Optimizer

• Replaces Mean-Variance optimization
• Minimizes Expected Shortfall (CVaR @ 95%)
• Uses simulated multi-path trajectories

---

💻 Proof of Concept (PoC)

1. EVT Residual Smoothing

import numpy as np
from scipy import stats

def apply_evt_smoothing(residuals, threshold_q=0.95):
    """
    Stabilizes residuals by replacing extreme outliers with GPD samples.
    """
    threshold = np.quantile(np.abs(residuals), threshold_q)
    extreme_mask = np.abs(residuals) > threshold
    
    if np.any(extreme_mask):
        exceedances = np.abs(residuals[extreme_mask]) - threshold
        shape, loc, scale = stats.genpareto.fit(exceedances)
        gpd_samples = stats.genpareto.rvs(shape, loc=loc, scale=scale, size=np.sum(extreme_mask))
        
        smoothed = residuals.copy()
        smoothed[extreme_mask] = np.sign(residuals[extreme_mask]) * (threshold + gpd_samples)
        return smoothed
    return residuals

---

2. GARCH-LSTM Joint Forecasting

import torch
import torch.nn as nn
from arch import arch_model
from sklearn.preprocessing import StandardScaler

class LSTMMean(nn.Module):
    def __init__(self, input_size=1, hidden_size=50):
        super().__init__()
        self.lstm = nn.LSTM(input_size, hidden_size, batch_first=True)
        self.fc = nn.Linear(hidden_size, 1)
    
    def forward(self, x):
        out, _ = self.lstm(x)
        return self.fc(out[:, -1, :])

def simulate_garch_lstm_joint(returns, corr_matrix, lookback, horizon, n_paths, device="cpu"):
    """
    Vectorized simulation using GARCH volatility, EVT residuals, and LSTM mean.
    """
    n_assets = returns.shape[1]
    garch_results = []
    lstm_models = []
    scalers = []
    all_stable_residuals = []
    
    # Phase 1: Training + EVT
    for i in range(n_assets):
        series = returns.iloc[:, i]
        res = arch_model(series * 100, vol="Garch", p=1, q=1).fit(disp="off")
        garch_results.append(res)
        smoothed_z = apply_evt_smoothing(res.std_resid)
        all_stable_residuals.append(np.sort(smoothed_z))
        
        scaler = StandardScaler().fit(series.values.reshape(-1, 1))
        scalers.append(scaler)
        lstm_models.append(LSTMMean().to(device))

    # Phase 2: Copula Sampling
    L = np.linalg.cholesky(corr_matrix)
    shocks = (np.random.normal(size=(horizon, n_paths, n_assets)) @ L.T)
    u_shocks = stats.norm.cdf(shocks)

    sim_returns = np.zeros((n_paths, horizon, n_assets))

    # Phase 3: Step-forward simulation (omitted for brevity)

    return sim_returns

---

3. CVaR Optimization

import cvxpy as cp

def optimize_cvar(simulated_returns, alpha=0.95, max_weight=0.3):
    """
    Minimizes Conditional Value at Risk (CVaR).
    """
    n_paths, n_assets = simulated_returns.shape
    weights = cp.Variable(n_assets)
    z = cp.Variable(n_paths)
    gamma = cp.Variable()
    
    loss = -simulated_returns @ weights
    objective = cp.Minimize(gamma + (1.0 / (n_paths * (1 - alpha))) * cp.sum(z))
    
    constraints = [
        z >= 0,
        z >= loss - gamma,
        cp.sum(weights) == 1,
        weights >= 0,
        weights <= max_weight
    ]
    
    prob = cp.Problem(objective, constraints)
    prob.solve(solver=cp.ECOS)
    
    return weights.value

I have done some sort of work on it which is curently working like a fine wine .If you want i can build a Qlib style wrapper that will have multiple try and except to reduce the error and also add different types of stochastic differential equation models in the wrapper that the user will just import and tune on the data

GARCH_LSTM_SVJ.py