옵션·선물 파생상품 기초: 그릭스 분석과 Python 헤징 시뮬레이션 실전 가이드

import numpy as np
from scipy.stats import norm
from dataclasses import dataclass


@dataclass
class BSMParams:
    """Black-Scholes-Merton 모델 파라미터"""
    S: float      # 기초자산 현재 가격
    K: float      # 행사가격
    T: float      # 만기까지 기간 (연 단위)
    r: float      # 무위험 이자율 (연율)
    sigma: float  # 변동성 (연율)
    q: float = 0  # 배당수익률


class BlackScholesModel:
    """Black-Scholes-Merton 옵션 가격 결정 모델"""

    def __init__(self, params: BSMParams):
        self.p = params
        self._d1 = self._calc_d1()
        self._d2 = self._calc_d2()

    def _calc_d1(self) -> float:
        p = self.p
        numerator = (
            np.log(p.S / p.K)
            + (p.r - p.q + 0.5 * p.sigma**2) * p.T
        )
        denominator = p.sigma * np.sqrt(p.T)
        return numerator / denominator

    def _calc_d2(self) -> float:
        return self._d1 - self.p.sigma * np.sqrt(self.p.T)

    def call_price(self) -> float:
        """유럽형 콜 옵션 가격"""
        p = self.p
        return (
            p.S * np.exp(-p.q * p.T) * norm.cdf(self._d1)
            - p.K * np.exp(-p.r * p.T) * norm.cdf(self._d2)
        )

    def put_price(self) -> float:
        """유럽형 풋 옵션 가격"""
        p = self.p
        return (
            p.K * np.exp(-p.r * p.T) * norm.cdf(-self._d2)
            - p.S * np.exp(-p.q * p.T) * norm.cdf(-self._d1)
        )

    def put_call_parity_check(self) -> dict:
        """풋-콜 패리티 검증: C - P = S*e^(-qT) - K*e^(-rT)"""
        p = self.p
        call = self.call_price()
        put = self.put_price()
        lhs = call - put
        rhs = p.S * np.exp(-p.q * p.T) - p.K * np.exp(-p.r * p.T)
        return {
            "call_price": round(call, 4),
            "put_price": round(put, 4),
            "C - P": round(lhs, 4),
            "S*e^(-qT) - K*e^(-rT)": round(rhs, 4),
            "parity_holds": abs(lhs - rhs) < 1e-10,
        }


# 사용 예시
params = BSMParams(S=100, K=105, T=0.5, r=0.05, sigma=0.2)
bsm = BlackScholesModel(params)
print(f"콜 옵션 가격: {bsm.call_price():.4f}")
print(f"풋 옵션 가격: {bsm.put_price():.4f}")
print(f"풋-콜 패리티: {bsm.put_call_parity_check()}")

import numpy as np
from scipy.stats import norm


class GreeksCalculator:
    """옵션 그릭스(1차, 2차) 계산기"""

    def __init__(self, S: float, K: float, T: float, r: float, sigma: float):
        self.S, self.K, self.T, self.r, self.sigma = S, K, T, r, sigma
        self.d1 = (
            (np.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
        )
        self.d2 = self.d1 - sigma * np.sqrt(T)

    def delta(self, option_type: str = "call") -> float:
        """Delta: 기초자산 가격 변화에 대한 옵션 가격의 민감도"""
        if option_type == "call":
            return norm.cdf(self.d1)
        return norm.cdf(self.d1) - 1

    def gamma(self) -> float:
        """Gamma: Delta의 변화율"""
        return norm.pdf(self.d1) / (self.S * self.sigma * np.sqrt(self.T))

    def theta(self, option_type: str = "call") -> float:
        """Theta: 시간 경과에 따른 옵션 가격 감소율 (1일 단위)"""
        common = -(self.S * norm.pdf(self.d1) * self.sigma) / (2 * np.sqrt(self.T))
        if option_type == "call":
            val = common - self.r * self.K * np.exp(-self.r * self.T) * norm.cdf(self.d2)
        else:
            val = common + self.r * self.K * np.exp(-self.r * self.T) * norm.cdf(-self.d2)
        return val / 365

    def vega(self) -> float:
        """Vega: 변동성 1% 변화에 대한 옵션 가격의 민감도"""
        return self.S * norm.pdf(self.d1) * np.sqrt(self.T) / 100

    def rho(self, option_type: str = "call") -> float:
        """Rho: 이자율 1% 변화에 대한 옵션 가격의 민감도"""
        if option_type == "call":
            return self.K * self.T * np.exp(-self.r * self.T) * norm.cdf(self.d2) / 100
        return -self.K * self.T * np.exp(-self.r * self.T) * norm.cdf(-self.d2) / 100

    def summary(self, option_type: str = "call") -> dict:
        return {
            "Delta": round(self.delta(option_type), 6),
            "Gamma": round(self.gamma(), 6),
            "Theta (daily)": round(self.theta(option_type), 6),
            "Vega (per 1%)": round(self.vega(), 6),
            "Rho (per 1%)": round(self.rho(option_type), 6),
        }


greeks = GreeksCalculator(S=100, K=100, T=0.25, r=0.05, sigma=0.2)
print("=== ATM 콜 옵션 그릭스 ===")
for name, value in greeks.summary("call").items():
    print(f"  {name}: {value}")

import numpy as np
import matplotlib.pyplot as plt


def plot_strategy_payoffs():
    """주요 헤징 전략의 만기 시 손익 다이어그램"""
    S = np.linspace(70, 130, 200)
    K, S0, C, P = 100, 100, 5, 4
    K_put, K_call = 95, 110

    fig, axes = plt.subplots(2, 2, figsize=(14, 10))

    # 1. Covered Call
    ax = axes[0, 0]
    stock_pnl = S - S0
    call_sold = np.where(S > K, -(S - K), 0) + C
    ax.plot(S, stock_pnl, '--', alpha=0.5, label='Stock Only')
    ax.plot(S, stock_pnl + call_sold, 'b-', linewidth=2, label='Covered Call')
    ax.axhline(y=0, color='k', linewidth=0.5)
    ax.set_title('Covered Call'); ax.legend(); ax.grid(True, alpha=0.3)

    # 2. Protective Put
    ax = axes[0, 1]
    put_bought = np.where(S < K, K - S, 0) - P
    ax.plot(S, stock_pnl, '--', alpha=0.5, label='Stock Only')
    ax.plot(S, stock_pnl + put_bought, 'g-', linewidth=2, label='Protective Put')
    ax.axhline(y=0, color='k', linewidth=0.5)
    ax.set_title('Protective Put'); ax.legend(); ax.grid(True, alpha=0.3)

    # 3. Collar
    ax = axes[1, 0]
    collar_put = np.where(S < K_put, K_put - S, 0) - 3
    collar_call = np.where(S > K_call, -(S - K_call), 0) + 2
    ax.plot(S, stock_pnl, '--', alpha=0.5, label='Stock Only')
    ax.plot(S, stock_pnl + collar_put + collar_call, 'r-', linewidth=2, label='Collar')
    ax.axhline(y=0, color='k', linewidth=0.5)
    ax.set_title('Collar Strategy'); ax.legend(); ax.grid(True, alpha=0.3)

    # 4. Straddle
    ax = axes[1, 1]
    straddle = np.maximum(S - K, 0) - C + np.maximum(K - S, 0) - P
    ax.plot(S, straddle, 'm-', linewidth=2, label='Long Straddle')
    ax.axhline(y=0, color='k', linewidth=0.5)
    ax.set_title('Long Straddle'); ax.legend(); ax.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig('hedging_payoff.png', dpi=150)
    plt.show()

plot_strategy_payoffs()

import numpy as np
import time


class MonteCarloOptionPricer:
    """Monte Carlo 시뮬레이션 기반 옵션 가격 결정"""

    def __init__(self, S0, K, T, r, sigma, n_sims=100_000, n_steps=252):
        self.S0, self.K, self.T, self.r, self.sigma = S0, K, T, r, sigma
        self.n_sims, self.n_steps = n_sims, n_steps

    def generate_paths(self) -> np.ndarray:
        """기하 브라운 운동(GBM)으로 주가 경로 생성 (안티테틱 변량법)"""
        dt = self.T / self.n_steps
        drift = (self.r - 0.5 * self.sigma**2) * dt
        diffusion = self.sigma * np.sqrt(dt)
        n_half = self.n_sims // 2
        Z = np.random.standard_normal((n_half, self.n_steps))
        Z = np.vstack([Z, -Z])  # Antithetic variates
        log_returns = drift + diffusion * Z
        log_paths = np.cumsum(log_returns, axis=1)
        log_paths = np.insert(log_paths, 0, 0, axis=1)
        return self.S0 * np.exp(log_paths)

    def european_price(self, option_type="call") -> dict:
        """유럽형 옵션 가격"""
        start = time.time()
        paths = self.generate_paths()
        S_T = paths[:, -1]
        if option_type == "call":
            payoffs = np.maximum(S_T - self.K, 0)
        else:
            payoffs = np.maximum(self.K - S_T, 0)
        df = np.exp(-self.r * self.T)
        price = df * np.mean(payoffs)
        se = df * np.std(payoffs) / np.sqrt(len(payoffs))
        return {
            "price": round(price, 4), "std_error": round(se, 4),
            "95_ci": [round(price - 1.96*se, 4), round(price + 1.96*se, 4)],
            "elapsed": round(time.time() - start, 3),
        }

    def asian_price(self, option_type="call") -> dict:
        """아시안 옵션 가격 (경로 의존)"""
        paths = self.generate_paths()
        avg_price = np.mean(paths[:, 1:], axis=1)
        if option_type == "call":
            payoffs = np.maximum(avg_price - self.K, 0)
        else:
            payoffs = np.maximum(self.K - avg_price, 0)
        df = np.exp(-self.r * self.T)
        price = df * np.mean(payoffs)
        se = df * np.std(payoffs) / np.sqrt(len(payoffs))
        return {"price": round(price, 4), "std_error": round(se, 4)}

    def barrier_price(self, barrier, barrier_type="down-and-out", option_type="call"):
        """배리어 옵션 가격"""
        paths = self.generate_paths()
        S_T = paths[:, -1]
        if barrier_type == "down-and-out":
            knocked = np.any(paths[:, 1:] <= barrier, axis=1)
        else:
            knocked = np.any(paths[:, 1:] >= barrier, axis=1)
        payoffs = np.maximum(S_T - self.K, 0) if option_type == "call" else np.maximum(self.K - S_T, 0)
        payoffs[knocked] = 0
        df = np.exp(-self.r * self.T)
        return {"price": round(df * np.mean(payoffs), 4), "knockout_ratio": round(np.mean(knocked), 4)}


# 실행 및 BSM과 비교
np.random.seed(42)
mc = MonteCarloOptionPricer(S0=100, K=100, T=0.5, r=0.05, sigma=0.2)
print("=== MC 유럽형 콜 ===", mc.european_price("call"))
print("=== 아시안 콜 ===", mc.asian_price("call"))
print("=== 배리어 콜 (Down-Out, B=85) ===", mc.barrier_price(85))

import numpy as np
from scipy.stats import norm


class DeltaHedgingSimulator:
    """Delta 헤징 포트폴리오의 성과를 시뮬레이션"""

    def __init__(self, S0, K, T, r, sigma, n_steps=63):
        self.S0, self.K, self.T, self.r, self.sigma = S0, K, T, r, sigma
        self.n_steps = n_steps
        self.dt = T / n_steps

    def _bsm_delta(self, S, tau):
        if tau <= 0:
            return 1.0 if S > self.K else 0.0
        d1 = (np.log(S/self.K) + (self.r + 0.5*self.sigma**2)*tau) / (self.sigma*np.sqrt(tau))
        return norm.cdf(d1)

    def _bsm_price(self, S, tau):
        if tau <= 0:
            return max(S - self.K, 0)
        d1 = (np.log(S/self.K) + (self.r + 0.5*self.sigma**2)*tau) / (self.sigma*np.sqrt(tau))
        d2 = d1 - self.sigma*np.sqrt(tau)
        return S*norm.cdf(d1) - self.K*np.exp(-self.r*tau)*norm.cdf(d2)

    def simulate(self, actual_sigma=None, rebalance_freq=1, seed=42):
        np.random.seed(seed)
        if actual_sigma is None:
            actual_sigma = self.sigma
        # 주가 경로 생성
        Z = np.random.standard_normal(self.n_steps)
        S = np.zeros(self.n_steps + 1)
        S[0] = self.S0
        for i in range(self.n_steps):
            S[i+1] = S[i] * np.exp(
                (self.r - 0.5*actual_sigma**2)*self.dt + actual_sigma*np.sqrt(self.dt)*Z[i]
            )
        # 헤징 실행
        premium = self._bsm_price(S[0], self.T)
        cash, shares = premium, 0.0
        for i in range(self.n_steps):
            tau = self.T - i * self.dt
            delta = self._bsm_delta(S[i], tau)
            if i % rebalance_freq == 0:
                trade = delta - shares
                cash -= trade * S[i]
                cash *= np.exp(self.r * self.dt * rebalance_freq)
                shares = delta
            else:
                cash *= np.exp(self.r * self.dt)
        # 만기 정산
        payoff = max(S[-1] - self.K, 0)
        portfolio = cash + shares * S[-1]
        pnl = portfolio - payoff
        return {
            "premium": round(premium, 4), "final_S": round(S[-1], 2),
            "payoff": round(payoff, 2), "hedge_pnl": round(pnl, 4),
            "error_pct": round(abs(pnl)/premium*100, 2),
        }


sim = DeltaHedgingSimulator(S0=100, K=100, T=0.25, r=0.05, sigma=0.2)

# 시나리오 비교
r1 = sim.simulate(actual_sigma=0.2, rebalance_freq=1)
r2 = sim.simulate(actual_sigma=0.3, rebalance_freq=1)
r3 = sim.simulate(actual_sigma=0.2, rebalance_freq=5)

print("=== 시나리오 1: 실현변동성=모델변동성 ===")
print(f"  프리미엄: {r1['premium']}, 헤징P&L: {r1['hedge_pnl']}, 오차: {r1['error_pct']}%")
print("=== 시나리오 2: 실현변동성(30%) > 모델변동성(20%) ===")
print(f"  프리미엄: {r2['premium']}, 헤징P&L: {r2['hedge_pnl']}, 오차: {r2['error_pct']}%")
print("=== 시나리오 3: 리밸런싱 주기 5일 ===")
print(f"  프리미엄: {r3['premium']}, 헤징P&L: {r3['hedge_pnl']}, 오차: {r3['error_pct']}%")