Risk Measurement and Stress Testing

import numpy as np
import pandas as pd

def historical_var(returns: pd.Series, confidence: float = 0.95, horizon: int = 1) -> float:
    """
    Args:
        returns: Daily return series
        confidence: Confidence level, commonly 0.95 or 0.99
        horizon: Holding period in days, default 1
    Returns:
        VaR value (positive means loss)
    """
    sorted_returns = returns.sort_values()
    index = int((1 - confidence) * len(sorted_returns))
    var_1d = -sorted_returns.iloc[index]
    return var_1d * np.sqrt(horizon)  # square-root-of-time rule

from scipy.stats import norm

def parametric_var(returns: pd.Series, confidence: float = 0.95, horizon: int = 1) -> float:
    mu = returns.mean()
    sigma = returns.std()
    z = norm.ppf(1 - confidence)
    var_1d = -(mu + z * sigma)
    return var_1d * np.sqrt(horizon)

def historical_cvar(returns: pd.Series, confidence: float = 0.95) -> float:
    """CVaR = the mean of all losses beyond VaR."""
    var = historical_var(returns, confidence)
    tail_losses = returns[returns < -var]
    return -tail_losses.mean() if len(tail_losses) > 0 else var

def max_drawdown_analysis(equity: pd.Series) -> dict:
    """
    Args:
        equity: Net-value series
    Returns:
        dict: max_drawdown, peak_date, trough_date, recovery_date, duration
    """
    peak = equity.cummax()
    drawdown = (equity - peak) / peak
    max_dd = drawdown.min()
    trough_idx = drawdown.idxmin()
    peak_idx = equity[:trough_idx].idxmax()

    # Recovery date
    recovery = equity[trough_idx:][equity[trough_idx:] >= equity[peak_idx]]
    recovery_date = recovery.index[0] if len(recovery) > 0 else None

    return {
        'max_drawdown': max_dd,
        'peak_date': peak_idx,
        'trough_date': trough_idx,
        'recovery_date': recovery_date,
        'underwater_days': (trough_idx - peak_idx).days,
        'recovery_days': (recovery_date - trough_idx).days if recovery_date else None
    }

def monte_carlo_gbm(S0: float, mu: float, sigma: float,
                     T: int = 252, n_paths: int = 10000) -> np.ndarray:
    """
    Args:
        S0: Initial price
        mu: Annualized return
        sigma: Annualized volatility
        T: Number of simulation days
        n_paths: Number of paths
    Returns:
        Price matrix of shape (n_paths, T)
    """
    dt = 1 / 252
    Z = np.random.standard_normal((n_paths, T))
    log_returns = (mu - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * Z
    prices = S0 * np.exp(np.cumsum(log_returns, axis=1))
    return prices

def analyze_mc_results(paths: np.ndarray, confidence: float = 0.95) -> dict:
    final_prices = paths[:, -1]
    returns = final_prices / paths[:, 0] - 1

    return {
        'mean_return': np.mean(returns),
        'median_return': np.median(returns),
        'std_return': np.std(returns),
        'var': -np.percentile(returns, (1 - confidence) * 100),
        'cvar': -np.mean(returns[returns < -np.percentile(returns, (1-confidence)*100)]),
        'prob_loss': np.mean(returns < 0),
        'worst_5pct': np.percentile(returns, 5),
        'best_5pct': np.percentile(returns, 95),
    }

STRESS_SCENARIOS = {
    'rate_shock_up_100bp': {
        'equity': -0.10,    # equities down 10%
        'bond_10y': -0.08,  # 10-year bonds down 8%
        'bond_2y': -0.02,   # short bonds down 2%
        'gold': +0.05,      # gold up 5%
        'btc': -0.15,       # BTC down 15%
    },
    'credit_crisis': {
        'equity': -0.25,
        'bond_10y': +0.05,  # government bonds act as a safe haven
        'credit_bond': -0.15,
        'gold': +0.10,
        'btc': -0.30,
    },
    'liquidity_dry_up': {
        'equity': -0.20,
        'bond_10y': -0.05,  # when liquidity is poor, everything falls
        'gold': -0.05,
        'btc': -0.40,
        'cash': 0.0,
    },
    'geopolitical_conflict': {
        'equity': -0.15,
        'bond_10y': +0.03,
        'gold': +0.15,
        'oil': +0.30,
        'btc': -0.20,
    },
}

from scipy.stats import genpareto

def fit_gpd_tail(returns: pd.Series, threshold_pct: float = 5.0) -> dict:
    """
    Fit the tail with a generalized Pareto distribution.
    Args:
        returns: Daily returns
        threshold_pct: Threshold percentile (take the worst X%)
    """
    threshold = np.percentile(returns, threshold_pct)
    exceedances = threshold - returns[returns < threshold]  # make positive

    # Fit GPD
    shape, loc, scale = genpareto.fit(exceedances)

    return {
        'threshold': threshold,
        'n_exceedances': len(exceedances),
        'shape_xi': shape,      # ξ>0 fat tail, ξ=0 exponential tail, ξ<0 bounded tail
        'scale_sigma': scale,
        'tail_type': 'fat tail (dangerous)' if shape > 0 else 'thin tail (safer)',
    }

Required:
- Return series (daily or higher frequency) or net-value series
- Portfolio weights (if it is a portfolio)

Optional:
- Benchmark returns (for relative risk analysis)
- Risk budget / constraint settings

## Risk Analysis Report

### Core Risk Metrics
| Metric | Value |
|------|-----|
| Daily volatility | 1.85% |
| Annualized volatility | 29.3% |
| Maximum drawdown | -32.5% (2024.09.15 → 2024.11.20) |
| VaR(95%, 1D) | -2.8% |
| CVaR(95%, 1D) | -4.2% |
| Skewness | -0.45 |
| Kurtosis | 5.2 (fat tail) |

### Stress-Test Results
| Scenario | Portfolio Loss | Stop Triggered |
|------|---------|----------|
| 2020 COVID replay | -18.5% | No |
| Rates +100bp | -12.3% | No |
| Liquidity dry-up | -28.7% | Yes |

### Monte Carlo Simulation (252 days, 10000 paths)
| Statistic | Value |
|------|-----|
| Expected return | +8.2% |
| Loss probability | 35% |
| Worst 5% scenario | -22.4% |

### Risk-Control Recommendations
1. Recommend setting a portfolio stop-loss at -15%
2. Tail risk is elevated; consider allocating 5% to gold as a hedge
3. Correlations rise in stressed markets, so diversification benefits will be discounted