ポートフォリオリスク管理：Pythonで実装するVaRとモンテカルロシミュレーション実践ガイド

VaR_α = -inf{x : P(L > x) <= 1 - α}

import numpy as np
import pandas as pd
from scipy import stats
import yfinance as yf
from datetime import datetime, timedelta

# ポートフォリオ銘柄とウェイト設定
tickers = ['AAPL', 'MSFT', 'GOOGL', 'AMZN', 'NVDA']
weights = np.array([0.25, 0.20, 0.20, 0.20, 0.15])
portfolio_value = 1_000_000_000  # 10億ウォン

# 株価データ収集（2年分）
end_date = datetime(2026, 3, 1)
start_date = end_date - timedelta(days=730)
prices = yf.download(tickers, start=start_date, end=end_date)['Adj Close']

# 日次対数収益率の計算
log_returns = np.log(prices / prices.shift(1)).dropna()

# ポートフォリオ収益率
portfolio_returns = log_returns.dot(weights)

# パラメトリックVaR計算
confidence_levels = [0.90, 0.95, 0.99]
mean_return = portfolio_returns.mean()
std_return = portfolio_returns.std()

print("=" * 60)
print("パラメトリックVaR（分散共分散法）")
print("=" * 60)
print(f"日次平均収益率: {mean_return:.6f}")
print(f"日次標準偏差:    {std_return:.6f}")
print(f"ポートフォリオ価値:  {portfolio_value:,.0f} ウォン")
print("-" * 60)

for cl in confidence_levels:
    z_score = stats.norm.ppf(1 - cl)
    var_pct = -(mean_return + z_score * std_return)
    var_amount = var_pct * portfolio_value
    print(f"{cl*100:.0f}% VaR: {var_pct:.4%} = {var_amount:,.0f} ウォン")

# 10日VaR（バーゼル規制基準）: sqrt(10)スケーリング
print("\n10日VaR（sqrt(T)スケーリング）:")
for cl in confidence_levels:
    z_score = stats.norm.ppf(1 - cl)
    var_1d = -(mean_return + z_score * std_return)
    var_10d = var_1d * np.sqrt(10)
    var_10d_amount = var_10d * portfolio_value
    print(f"{cl*100:.0f}% 10日VaR: {var_10d:.4%} = {var_10d_amount:,.0f} ウォン")

# 共分散行列ベースのポートフォリオVaR
cov_matrix = log_returns.cov()

# ポートフォリオ分散 = w^T * Σ * w
portfolio_variance = weights @ cov_matrix.values @ weights
portfolio_std = np.sqrt(portfolio_variance)

# 個別資産のVaR寄与度（Component VaR）
marginal_var = cov_matrix.values @ weights / portfolio_std
component_var = weights * marginal_var

print("\n" + "=" * 60)
print("コンポーネントVaR分析（99%信頼水準）")
print("=" * 60)
z_99 = stats.norm.ppf(0.99)

for i, ticker in enumerate(tickers):
    cvar_pct = component_var[i] * z_99
    cvar_amount = cvar_pct * portfolio_value
    contribution_pct = component_var[i] / portfolio_std * 100
    print(f"{ticker}: ウェイト {weights[i]*100:.0f}% | "
          f"リスク寄与度 {contribution_pct:.1f}% | "
          f"Component VaR {cvar_amount:,.0f} ウォン")

# 分散効果の確認
individual_vars = np.array([
    weights[i] * np.sqrt(cov_matrix.values[i, i]) * z_99
    for i in range(len(tickers))
])
undiversified_var = individual_vars.sum() * portfolio_value
diversified_var = portfolio_std * z_99 * portfolio_value
diversification_benefit = undiversified_var - diversified_var

print(f"\n非分散VaR（合算）:   {undiversified_var:,.0f} ウォン")
print(f"分散VaR（ポートフォリオ）: {diversified_var:,.0f} ウォン")
print(f"分散効果:            {diversification_benefit:,.0f} ウォン "
      f"({diversification_benefit/undiversified_var*100:.1f}% 削減)")

def historical_var(returns, confidence_level=0.95, portfolio_value=1_000_000_000):
    """
    ヒストリカルシミュレーション方式でVaRを計算する。
    """
    sorted_returns = returns.sort_values()
    n = len(sorted_returns)

    percentile_index = int(np.floor((1 - confidence_level) * n))
    var_pct = -sorted_returns.iloc[percentile_index]

    var_pct_interp = -np.percentile(returns, (1 - confidence_level) * 100)

    return {
        'var_pct': var_pct,
        'var_pct_interp': var_pct_interp,
        'var_amount': var_pct * portfolio_value,
        'var_amount_interp': var_pct_interp * portfolio_value,
        'n_observations': n,
        'worst_loss': -sorted_returns.iloc[0],
        'worst_loss_date': sorted_returns.index[0],
    }

results = {}
for cl in [0.90, 0.95, 0.99]:
    result = historical_var(portfolio_returns, cl)
    results[cl] = result
    print(f"\n{cl*100:.0f}% ヒストリカルVaR:")
    print(f"  パーセンタイル方式: {result['var_pct']:.4%} = {result['var_amount']:,.0f} ウォン")
    print(f"  補間法方式:   {result['var_pct_interp']:.4%} = {result['var_amount_interp']:,.0f} ウォン")

print(f"\n観測期間内最大損失: {results[0.95]['worst_loss']:.4%} "
      f"({results[0.95]['worst_loss_date'].strftime('%Y-%m-%d')})")

rolling_window = 250
rolling_var_95 = portfolio_returns.rolling(rolling_window).apply(
    lambda x: -np.percentile(x, 5), raw=True
)

print(f"\nローリングVaR（250日ウィンドウ）統計:")
print(f"  平均 95% VaR: {rolling_var_95.mean():.4%}")
print(f"  最大 95% VaR: {rolling_var_95.max():.4%}（最もリスクが高い時期）")
print(f"  最小 95% VaR: {rolling_var_95.min():.4%}（最もリスクが低い時期）")

def monte_carlo_var(
    returns_df,
    weights,
    portfolio_value,
    n_simulations=50_000,
    n_days=1,
    confidence_level=0.95,
    seed=42
):
    """
    モンテカルロシミュレーションでポートフォリオVaRを計算する。

    幾何ブラウン運動（GBM）ベースで相関関係が反映された
    多変量正規分布からシナリオを生成する。
    """
    np.random.seed(seed)

    mean_returns = returns_df.mean().values
    cov_matrix = returns_df.cov().values
    n_assets = len(weights)

    L = np.linalg.cholesky(cov_matrix)

    portfolio_sim_returns = np.zeros(n_simulations)

    for sim in range(n_simulations):
        cumulative_return = 0
        for day in range(n_days):
            Z = np.random.standard_normal(n_assets)
            correlated_returns = mean_returns + L @ Z
            daily_port_return = weights @ correlated_returns
            cumulative_return += daily_port_return

        portfolio_sim_returns[sim] = cumulative_return

    sim_pnl = portfolio_sim_returns * portfolio_value

    var_pct = -np.percentile(portfolio_sim_returns, (1 - confidence_level) * 100)
    var_amount = var_pct * portfolio_value

    var_threshold = np.percentile(portfolio_sim_returns, (1 - confidence_level) * 100)
    tail_losses = portfolio_sim_returns[portfolio_sim_returns <= var_threshold]
    cvar_pct = -tail_losses.mean()
    cvar_amount = cvar_pct * portfolio_value

    return {
        'var_pct': var_pct,
        'var_amount': var_amount,
        'cvar_pct': cvar_pct,
        'cvar_amount': cvar_amount,
        'sim_returns': portfolio_sim_returns,
        'sim_pnl': sim_pnl,
        'n_simulations': n_simulations,
        'n_days': n_days,
        'confidence_level': confidence_level,
    }

mc_results = {}
for cl in [0.90, 0.95, 0.99]:
    result = monte_carlo_var(
        log_returns, weights, portfolio_value,
        n_simulations=100_000,
        n_days=1,
        confidence_level=cl,
        seed=42
    )
    mc_results[cl] = result
    print(f"\n{cl*100:.0f}% モンテカルロVaR（100,000シミュレーション）:")
    print(f"  VaR:  {result['var_pct']:.4%} = {result['var_amount']:,.0f} ウォン")
    print(f"  CVaR: {result['cvar_pct']:.4%} = {result['cvar_amount']:,.0f} ウォン")

print("\nシミュレーション数による95% VaRの収束:")
for n_sim in [1_000, 5_000, 10_000, 50_000, 100_000, 500_000]:
    result = monte_carlo_var(
        log_returns, weights, portfolio_value,
        n_simulations=n_sim, confidence_level=0.95, seed=42
    )
    print(f"  N={n_sim:>7,}: VaR = {result['var_pct']:.6%}")

CVaR_α = E[L | L > VaR_α]

def comprehensive_risk_metrics(returns, portfolio_value, confidence_level=0.95):
    """
    ポートフォリオの総合リスク指標を計算する。

    VaR、CVaR以外にも最大ドローダウン（MDD）、ボラティリティ、
    歪度、尖度などを合わせて計算し、リスクの全体像を提供する。
    """
    sorted_returns = returns.sort_values()
    n = len(sorted_returns)

    var_pct = -np.percentile(returns, (1 - confidence_level) * 100)

    var_threshold = np.percentile(returns, (1 - confidence_level) * 100)
    tail_returns = returns[returns <= var_threshold]
    cvar_pct = -tail_returns.mean()

    daily_vol = returns.std()
    annual_vol = daily_vol * np.sqrt(252)

    skewness = returns.skew()
    kurtosis = returns.kurtosis()

    cumulative = (1 + returns).cumprod()
    running_max = cumulative.cummax()
    drawdown = (cumulative - running_max) / running_max
    max_drawdown = drawdown.min()

    parametric_var = -(returns.mean() + stats.norm.ppf(1 - confidence_level) * returns.std())
    tail_ratio = var_pct / parametric_var if parametric_var > 0 else float('inf')

    results = {
        'VaR': var_pct,
        'CVaR': cvar_pct,
        'VaR（金額）': var_pct * portfolio_value,
        'CVaR（金額）': cvar_pct * portfolio_value,
        'CVaR/VaR比率': cvar_pct / var_pct if var_pct > 0 else 0,
        '日次ボラティリティ': daily_vol,
        '年率ボラティリティ': annual_vol,
        '歪度': skewness,
        '超過尖度': kurtosis,
        '最大ドローダウン': max_drawdown,
        'テール比率（実際/パラメトリック）': tail_ratio,
    }

    print("=" * 60)
    print(f"総合リスク分析（{confidence_level*100:.0f}%信頼水準）")
    print("=" * 60)
    for key, value in results.items():
        if '金額' in key:
            print(f"  {key}: {value:,.0f} ウォン")
        elif '比率' in key or 'ボラティリティ' in key or 'ドローダウン' in key:
            print(f"  {key}: {value:.4%}")
        else:
            print(f"  {key}: {value:.6f}")

    print("\n[解釈ガイド]")
    if kurtosis > 1:
        print(f"  - 超過尖度 {kurtosis:.2f}: 正規分布よりテールが厚い"
              f"（ファットテール）。パラメトリックVaRがリスクを過小評価する可能性が高い。")
    if skewness < -0.5:
        print(f"  - 歪度 {skewness:.2f}: 負の非対称性。"
              f"大きな損失が大きな利益よりも頻繁に/大きく発生する傾向。")
    if tail_ratio > 1.2:
        print(f"  - テール比率 {tail_ratio:.2f}: 実際のテールリスクが"
              f"正規分布仮定対比 {(tail_ratio-1)*100:.0f}% 大きい。")

    return results

risk_metrics = comprehensive_risk_metrics(portfolio_returns, portfolio_value, 0.95)

from pypfopt import EfficientFrontier, risk_models, expected_returns
from pypfopt import HRPOpt
from pypfopt.efficient_frontier import EfficientCVaR

mu = expected_returns.mean_historical_return(prices)
S = risk_models.sample_cov(prices)

ef_minvol = EfficientFrontier(mu, S)
ef_minvol.min_volatility()
w_minvol = ef_minvol.clean_weights()
perf_minvol = ef_minvol.portfolio_performance(verbose=False)

ef_sharpe = EfficientFrontier(mu, S)
ef_sharpe.max_sharpe(risk_free_rate=0.035)
w_sharpe = ef_sharpe.clean_weights()
perf_sharpe = ef_sharpe.portfolio_performance(verbose=False)

hrp = HRPOpt(log_returns)
hrp.optimize()
w_hrp = hrp.clean_weights()
perf_hrp = hrp.portfolio_performance(verbose=False)

ef_cvar = EfficientCVaR(mu, log_returns)
ef_cvar.min_cvar()
w_cvar = ef_cvar.clean_weights()

print("=" * 70)
print("ポートフォリオ最適化結果比較")
print("=" * 70)
print(f"{'戦略':<20} {'期待収益率':>10} {'ボラティリティ':>10} {'シャープ':>8}")
print("-" * 70)
print(f"{'最小ボラティリティ':<20} {perf_minvol[0]:>10.2%} {perf_minvol[1]:>10.2%} "
      f"{perf_minvol[2]:>8.3f}")
print(f"{'最大シャープ':<20} {perf_sharpe[0]:>10.2%} {perf_sharpe[1]:>10.2%} "
      f"{perf_sharpe[2]:>8.3f}")
print(f"{'HRP':<20} {perf_hrp[0]:>10.2%} {perf_hrp[1]:>10.2%} "
      f"{perf_hrp[2]:>8.3f}")

print("\nウェイト比較:")
for ticker in tickers:
    print(f"  {ticker}: 最小ボラ {w_minvol.get(ticker,0):.1%} | "
          f"最大シャープ {w_sharpe.get(ticker,0):.1%} | "
          f"HRP {w_hrp.get(ticker,0):.1%} | "
          f"最小CVaR {w_cvar.get(ticker,0):.1%}")

def backtest_var(returns, var_series, confidence_level=0.95, method_name=""):
    """
    VaRバックテスティング：Kupiec検定の実施
    """
    breaches = returns < -var_series
    n_breaches = breaches.sum()
    n_total = len(returns)
    breach_rate = n_breaches / n_total
    expected_rate = 1 - confidence_level
    expected_breaches = expected_rate * n_total

    print(f"\n{'=' * 60}")
    print(f"VaRバックテスティング結果: {method_name}")
    print(f"{'=' * 60}")
    print(f"  観測期間:       {n_total} 営業日")
    print(f"  VaR超過回数:   {n_breaches} 回")
    print(f"  期待超過回数:  {expected_breaches:.1f} 回")
    print(f"  実際超過比率:  {breach_rate:.2%}")
    print(f"  期待超過比率:  {expected_rate:.2%}")

    p = expected_rate
    p_hat = breach_rate if breach_rate > 0 else 1e-10

    if n_breaches == 0:
        lr_stat = -2 * (n_total * np.log(1 - p) - n_total * np.log(1 - p_hat))
    elif n_breaches == n_total:
        lr_stat = -2 * (n_total * np.log(p) - n_total * np.log(p_hat))
    else:
        lr_stat = -2 * (
            n_breaches * np.log(p) +
            (n_total - n_breaches) * np.log(1 - p) -
            n_breaches * np.log(p_hat) -
            (n_total - n_breaches) * np.log(1 - p_hat)
        )

    p_value = 1 - stats.chi2.cdf(lr_stat, df=1)

    print(f"\n  Kupiec LR統計量: {lr_stat:.4f}")
    print(f"  p値:          {p_value:.4f}")

    if p_value > 0.05:
        verdict = "合格（モデル適合）"
    else:
        verdict = "不合格（モデル不適合 - 再検討必要）"
    print(f"  判定（5%有意水準）: {verdict}")

    if n_total >= 250:
        scale_factor = n_total / 250
        scaled_breaches = n_breaches / scale_factor
        if scaled_breaches <= 4:
            zone = "グリーン（Green Zone） - 良好"
        elif scaled_breaches <= 9:
            zone = "イエロー（Yellow Zone） - 注意必要"
        else:
            zone = "レッド（Red Zone） - モデル修正必須"
        print(f"  バーゼル信号灯:      {zone}")

    if n_breaches > 0:
        consecutive = 0
        max_consecutive = 0
        for b in breaches:
            if b:
                consecutive += 1
                max_consecutive = max(max_consecutive, consecutive)
            else:
                consecutive = 0
        print(f"\n  最大連続超過:   {max_consecutive} 営業日")
        if max_consecutive >= 3:
            print(f"  警告: 連続超過 {max_consecutive}回発生。"
                  f"ボラティリティクラスタリング未反映の可能性。")

    return {
        'n_breaches': n_breaches,
        'breach_rate': breach_rate,
        'lr_stat': lr_stat,
        'p_value': p_value,
        'pass': p_value > 0.05,
    }

window = 250
test_returns = portfolio_returns.iloc[window:]

rolling_mean = portfolio_returns.rolling(window).mean().iloc[window:]
rolling_std = portfolio_returns.rolling(window).std().iloc[window:]
z_95 = stats.norm.ppf(0.95)
parametric_var_series = -(rolling_mean - z_95 * rolling_std)

historical_var_series = portfolio_returns.rolling(window).apply(
    lambda x: -np.percentile(x, 5), raw=True
).iloc[window:]

bt_param = backtest_var(test_returns, parametric_var_series, 0.95, "パラメトリックVaR")
bt_hist = backtest_var(test_returns, historical_var_series, 0.95, "ヒストリカルVaR")