포트폴리오 리스크 관리: 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 * Sigma * 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를 계산한다.

    Parameters
    ----------
    returns : pd.Series
        포트폴리오 일간 수익률 시계열
    confidence_level : float
        신뢰수준 (예: 0.95, 0.99)
    portfolio_value : float
        포트폴리오 총 가치 (원)

    Returns
    -------
    dict
        VaR 계산 결과
    """
    sorted_returns = returns.sort_values()
    n = len(sorted_returns)

    # 하위 (1-confidence_level) 백분위수
    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],
    }

# 히스토리컬 VaR 계산
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')})")

# 히스토리컬 VaR의 롤링 분석
rolling_window = 250  # 1년 윈도우
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) 기반으로 상관관계가 반영된
    다변량 정규분포에서 시나리오를 생성한다.

    Parameters
    ----------
    returns_df : pd.DataFrame
        개별 자산의 일간 로그 수익률 DataFrame
    weights : np.ndarray
        자산 비중 벡터
    portfolio_value : float
        포트폴리오 총 가치 (원)
    n_simulations : int
        시뮬레이션 횟수
    n_days : int
        보유 기간 (일)
    confidence_level : float
        신뢰수준
    seed : int
        난수 시드 (재현성 보장)

    Returns
    -------
    dict
        VaR, CVaR, 시뮬레이션 결과 등
    """
    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 계산
    var_pct = -np.percentile(portfolio_sim_returns, (1 - confidence_level) * 100)
    var_amount = var_pct * portfolio_value

    # CVaR (Expected Shortfall) 계산
    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,
    }

# 몬테카를로 VaR 실행
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} 원")

# 시뮬레이션 수에 따른 VaR 수렴 확인
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
    var_pct = -np.percentile(returns, (1 - confidence_level) * 100)

    # CVaR (Expected Shortfall)
    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()  # 초과 첨도 (정규분포 = 0)

    # 최대 낙폭 (Maximum Drawdown)
    cumulative = (1 + returns).cumprod()
    running_max = cumulative.cummax()
    drawdown = (cumulative - running_max) / running_max
    max_drawdown = drawdown.min()

    # 정규분포 대비 VaR 비율 (테일 위험 지표)
    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)

# 1. 평균-분산 최적화 (최소 변동성 포트폴리오)
ef_minvol = EfficientFrontier(mu, S)
ef_minvol.min_volatility()
w_minvol = ef_minvol.clean_weights()
perf_minvol = ef_minvol.portfolio_performance(verbose=False)

# 2. 최대 샤프 비율 포트폴리오
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)

# 3. 계층적 리스크 패리티 (HRP)
hrp = HRPOpt(log_returns)
hrp.optimize()
w_hrp = hrp.clean_weights()
perf_hrp = hrp.portfolio_performance(verbose=False)

# 4. 최소 CVaR 포트폴리오
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 검정 수행

    Parameters
    ----------
    returns : pd.Series
        실제 포트폴리오 수익률
    var_series : pd.Series
        동일 기간의 VaR 추정값 (양수)
    confidence_level : float
        VaR 신뢰수준
    method_name : str
        방법론 이름 (출력용)
    """
    # VaR 초과 여부 판정
    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%}")

    # Kupiec 우도비 검정
    p = expected_rate
    p_hat = breach_rate if breach_rate > 0 else 1e-10

    # LR 통계량 계산
    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)
        )

    # 카이제곱 분포 (자유도 1)와 비교
    p_value = 1 - stats.chi2.cdf(lr_stat, df=1)

    print(f"\n  Kupiec LR 통계량: {lr_stat:.4f}")
    print(f"  p-value:          {p_value:.4f}")

    if p_value > 0.05:
        verdict = "통과 (모델 적합)"
    else:
        verdict = "실패 (모델 부적합 - 재검토 필요)"
    print(f"  판정 (5% 유의수준): {verdict}")

    # 바젤 신호등 체계 (Traffic Light)
    # 250일 기준: 초록(0~4), 노랑(5~9), 빨강(10+)
    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}")

    # 연속 초과 분석 (Christoffersen 독립성 예비 분석)
    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,
    }

# 롤링 VaR로 백테스팅 수행
window = 250
test_returns = portfolio_returns.iloc[window:]

# 파라메트릭 VaR 롤링 계산
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)

# 히스토리컬 VaR 롤링 계산
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")