VaR与CVaR实战：量化风险管理的双剑合璧

import numpy as np
import pandas as pd

def historical_var(returns, confidence=0.99):
    """
    历史模拟法计算VaR
    returns: 收益率序列
    confidence: 置信水平
    """
    var = np.percentile(returns, (1 - confidence) * 100)
    return abs(var)

# 示例：计算沪深300的99% VaR
hs300_returns = pd.read_csv('hs300_returns.csv')['return']
var_99 = historical_var(hs300_returns, 0.99)
print(f"99% VaR: {var_99:.2%}")

def monte_carlo_var(returns, confidence=0.99, num_simulations=10000):
    """
    蒙特卡洛模拟计算VaR
    使用GARCH(1,1)建模波动率聚集
    """
    from arch import arch_model
    
    # 拟合GARCH模型
    model = arch_model(returns * 100, vol='Garch', p=1, q=1)
    model_fit = model.fit(disp='off')
    
    # 模拟未来收益率
    forecasts = model_fit.forecast(horizon=1)
    simulated_returns = np.random.normal(
        loc=0,
        scale=forecasts.variance.values[-1, 0]**0.5 / 100,
        size=num_simulations
    )
    
    var = np.percentile(simulated_returns, (1 - confidence) * 100)
    return abs(var)

var_99_mc = monte_carlo_var(hs300_returns, 0.99)
print(f"99% VaR (蒙特卡洛): {var_99_mc:.2%}")

def multi_horizon_risk(returns, horizons=[1, 5, 10, 22]):
    """
    计算多个时间维度的VaR和CVaR
    horizons: 交易日数（1天、1周、2周、1个月）
    """
    results = []
    
    for days in horizons:
        # 计算滚动窗口的VaR和CVaR
        rolling_var = returns.rolling(window=252).apply(
            lambda x: historical_var(x, 0.99)
        )
        rolling_cvar = returns.rolling(window=252).apply(
            lambda x: historical_cvar(x, 0.99)
        )
        
        results.append({
            'horizon': f"{days}天",
            'VaR_99': rolling_var.iloc[-1],
            'CVaR_99': rolling_cvar.iloc[-1]
        })
    
    return pd.DataFrame(results)

risk_table = multi_horizon_risk(returns)
print(risk_table)

def backtest_var(returns, var_series, confidence=0.99):
    """
    回测VaR模型：检查实际突破次数是否符合预期
    """
    expected_breach_rate = 1 - confidence
    actual_breaches = (returns < -var_series).sum()
    breach_rate = actual_breaches / len(returns)
    
    print(f"期望突破率: {expected_breach_rate:.2%}")
    print(f"实际突破率: {breach_rate:.2%}")
    print(f"突破次数: {actual_breaches}")
    
    # 使用Kupiec检验（似然比检验）
    from scipy.stats import chi2
    lr_stat = -2 * np.log(
        (1 - expected_breach_rate)**(len(returns) - actual_breaches) *
        expected_breach_rate**actual_breaches /
        ((1 - breach_rate)**(len(returns) - actual_breaches) *
         breach_rate**actual_breaches)
    )
    p_value = 1 - chi2.cdf(lr_stat, df=1)
    print(f"Kupiec检验 p-value: {p_value:.4f}")
    
    if p_value < 0.05:
        print("❌ VaR模型拒绝原假设（模型不准确）")
    else:
        print("✅ VaR模型通过回测")

backtest_var(returns, rolling_var, 0.99)

import matplotlib.pyplot as plt
import matplotlib.dates as mdates

def plot_risk_dashboard(returns, var_series, cvar_series):
    fig, axes = plt.subplots(2, 1, figsize=(14, 10))
    
    # 图1：收益率 + VaR/CVaR带
    axes[0].plot(returns.index[-252:], returns[-252:], 
                label='日收益率', alpha=0.7)
    axes[0].fill_between(
        returns.index[-252:],
        -var_series[-252:], 0,
        alpha=0.3, color='orange',
        label='99% VaR'
    )
    axes[0].fill_between(
        returns.index[-252:],
        -cvar_series[-252:], -var_series[-252:],
        alpha=0.2, color='red',
        label='99% CVaR (尾部)'
    )
    axes[0].axhline(y=0, color='black', linestyle='--', alpha=0.5)
    axes[0].set_title('投资组合风险度量：VaR与CVaR')
    axes[0].legend()
    
    # 图2：突破次数累积图
    breaches = (returns < -var_series).astype(int)
    cumulative_breaches = breaches.cumsum()
    axes[1].plot(cumulative_breaches.index[-252:], 
                cumulative_breaches[-252:],
                label='VaR突破次数', color='red')
    axes[1].axhline(y=len(returns[-252:]) * 0.01, 
                   color='green', linestyle='--',
                   label='期望突破次数 (1%)')
    axes[1].set_title('VaR突破次数累积图')
    axes[1].legend()
    
    plt.tight_layout()
    plt.savefig('var_cvar_dashboard.png', dpi=300, bbox_inches='tight')
    plt.show()

plot_risk_dashboard(returns, rolling_var, rolling_cvar)

def historical_stress_test(portfolio_returns, stress_periods):
    """
    用历史事件重演测试投资组合
    stress_periods: [(开始日期, 结束日期, 事件名称), ...]
    """
    results = []
    
    for start, end, event in stress_periods:
        period_returns = portfolio_returns[start:end]
        cumulative_return = (1 + period_returns).prod() - 1
        max_drawdown = calculate_max_drawdown(period_returns)
        
        results.append({
            '事件': event,
            '期间收益率': f"{cumulative_return:.2%}",
            '最大回撤': f"{max_drawdown:.2%}"
        })
    
    return pd.DataFrame(results)

# 定义压力场景
stress_scenarios = [
    ('2015-06-01', '2015-08-31', '2015年股灾'),
    ('2018-01-01', '2018-12-31', '2018年贸易战'),
    ('2020-01-01', '2020-03-31', '2020年疫情崩盘'),
    ('2022-01-01', '2022-04-30', '2022年上海封城')
]

stress_results = historical_stress_test(returns, stress_scenarios)
print(stress_results)

def monte_carlo_stress_test(returns, num_simulations=10000, horizon=22):
    """
    蒙特卡洛压力测试：模拟未来1个月的极端情况
    """
    # 使用历史模拟法生成模拟收益率
    simulated_returns = np.random.choice(
        returns, 
        size=(num_simulations, horizon),
        replace=True
    )
    
    # 计算每轮模拟的累积收益率
    cumulative_returns = (1 + simulated_returns).prod(axis=1) - 1
    
    # 找出最坏的5%情况
    worst_5_percent = np.percentile(cumulative_returns, 5)
    
    print(f"未来1个月，有5%概率损失超过: {abs(worst_5_percent):.2%}")
    print(f"最坏情况下（0.1%分位数）: {abs(np.percentile(cumulative_returns, 0.1)):.2%}")
    
    # 可视化
    plt.figure(figsize=(10, 6))
    plt.hist(cumulative_returns, bins=50, alpha=0.7, edgecolor='black')
    plt.axvline(x=worst_5_percent, color='red', 
               linestyle='--', label='5%分位数')
    plt.xlabel('1个月累积收益率')
    plt.ylabel('频率')
    plt.title('蒙特卡洛压力测试结果（10000次模拟）')
    plt.legend()
    plt.show()

monte_carlo_stress_test(returns)