미래주가의 분포를 활용하여 콜옵션 델타 구하기

 

 

 

이번 글에서는 GBM 모델로 생성된 주가의 누적분포함수, 확률밀도함수

GBM 모델로 생성된 주가의 누적분포함수, 확률밀도함수

에서 다룬 내용을 활용하여,  파생상품의 민감도를 구하는 방법을 설명해 보겠습니다.

 

먼저 잠시 복습을 해볼까요?

 

 

복습

 

GBM 모델을 따르는 기초자산의 특정 만기 시점 $T$에서의 가격인 $S_T$의 누적분포함수(cdf)와 확률밀도 함수(pdf)는 아래와 같습니다.

 

기초자산  $S_t$가 GBM 모델 

$$dS_t/S_t= (r-q) dt +\sigma dW_t$$

을 따를 때($r$는 무위험 이자율, $q$는 연속배당률,  $\sigma$는 변동성(volatility), 그리고 $W_t$는 위너프로세스)

 

만일 현재시점의 주가를 $S$라 하면, $S_T$는

$$\ln(S_T/S) = (r-q-\frac12\sigma^2)T +\sigma W_T$$ 를 만족하고

 

$S_T$의	수식
누적분포함수 $F(x)$	$$F(x)=\Phi(d(x))$$ 여기서 $\Phi(\cdot)$는 표준정규분포의 누적분포함수이고 $$d(x) = \frac{\ln(x/S)-(r-q-\frac12\sigma^2)T}{\sigma\sqrt{T}}$$
확률밀도함수 $g(x)$	$$g(x)  = \frac{1}{x\sigma\sqrt{T}} \phi(d(x)),$$ 여기서 $\phi(\cdot)$는 표준정규분포의 확률밀도함수

 

입니다. 이 사실을 이용하여 유러피안 콜옵션(European Call)의 델타를 구해보기로 하죠.

 

 

 

콜옵션의 델타

 

만기 $T$이고 행사가 $K$인 유로피안 콜옵션의 가격 $c$는

$$c(S) = e^{-rT} \mathbb{E}(\max(S_T-K,0)) \tag{1}$$
입니다(현재시점 $t=0$에서의 옵션 가격으로 간주하여, $c$를 $S$에 대한 1 변수 함수로 표현하였습니다.)

위의 복습 파트에서, $S_T$의 확률밀도 함수를 $g(\cdot)$라 하였으므로 식(1)은

$$\begin{align}  c(S) & = e^{-rT} \int_0^\infty \max(x-K,0) g(x)dx\\       & =  \int_0^\infty e^{-rT}\max(x-K,0) g(x) \tag{2}dx \end{align}$$

즉 식(1)의 기댓값을 확률밀도함수를 이용하여 다시 쓴 것입니다. 변수 $x$는 $S_T$를 의미하고 $x$의 범위가 $(0,\infty)$이므로 적분 구간을 위와 같이 한 것입니다. 이제 위 콜옵션의 델타를 구해보겠습니다.

 

델타 구하기

 

식(2)을 $S$로 미분하여 델타를 구합니다.

$$\begin{align} \frac{\partial c}{\partial S} & = \frac{\partial}{\partial S} \int_0^\infty e^{-rT}\max(x-K,0) g(x) dx\\ & =  \int_0^\infty \frac{\partial}{\partial S} \left( e^{-rT}\max(x-K,0) g(x) \right)dx\\ & =  \int_0^\infty e^{-rT}\max(x-K,0) \frac{\partial g(x)}{\partial S} dx\\ & =  \int_0^\infty e^{-rT}\max(x-K,0) \frac{\frac{\partial g(x)}{\partial S}}{g(x)} \cdot g(x) dx\\ & = \mathbb{E} \left(e^{-rT}\max(x-K,0) \frac{\frac{\partial g(x)}{\partial S}}{g(x)} \right) \\ & = \mathbb{E} \left(e^{-rT}\max(x-K,0) \frac{\partial \ln(g(x))}{\partial S} \right)\tag{3} \end{align}$$

 

위에서 4번째 등식은 $g(x)$로 나눠주고, 다시 $g(x)$를 곱해 준 것입니다. 또  chain rule에 의해
$$\frac{\partial \ln(g(x))}{\partial S} = \frac{1}{g(x)} \cdot \frac{\partial g}{\partial S}$$ 가 성립하기 때문에 6번째 등식처럼 쓸 수 있습니다.

$$g(x) = \frac{1}{x\sigma\sqrt{T}} \phi\left( d(x)\right)~,~ d(x)=\frac{\ln(x/S)-(r-q-\frac12\sigma^2)T}{\sigma\sqrt{T}}$$
이므로

$$\begin{align} \frac{\partial}{\partial S} \ln g(x) & = \frac{\partial}{\partial S} \left[ \ln \left(\frac{1}{x\sigma\sqrt{T}}\phi(d(x)) \right)\right]\\ & = \frac{\partial}{\partial S} \left[ -\ln x -\ln (\sigma\sqrt{T}) +\ln \phi(d(x)) \right]\\ & = 0 + 0 + \frac{\partial}{\partial S} \left[ \ln \phi(d(x)) \right]\\ & = \frac{\phi'(d(x))}{\phi(d(x))} \frac{\partial d(x)}{\partial S}\\ & = -d(x) \cdot \left( - \frac{1}{S\sigma\sqrt{T}} \right)\\ & = \frac{d(x)}{S\sigma\sqrt{T}} \tag{4} \end{align}$$
( 위 식 5번째 등식은  $\phi'(z) =-z \phi(z)$의 사실이 쓰임)

따라서 식(3), (4)를 묶어서 쓰면

 

$$\frac{\partial c}{\partial S} = \mathbb{E} \left(e^{-rT}\max(x-K,0)\frac{d(x)}{S\sigma\sqrt{T}}  \right)\tag{5}$$

 

그렇다면 식(5)을 어떻게 구할 수 있을까요?

식(5) 우변의 $x$는 현재가가 $S$인 기초자산의 만기 시 가격 분포를 나타내는 확률변수입니다. 즉, $x=S_T$이죠. $S_T$는 어떤 표준정규분포난수 $z$에 대해

 

$$S_T = S \exp \left((r-q-\frac12\sigma^2)T+\sigma\sqrt{T}z \right)$$

의 형태로 생성이 되죠. 이를 다시 함수 $d(\cdot)$에 집어넣으면,

$$d(S_T) = \frac{\ln(S_T/S)-(r-q-\frac12\sigma^2)T}{\sigma\sqrt{T}} =z$$ 가 됩니다. 따라서!! 식 (5)는

 

$$\frac{\partial c}{\partial S} = \mathbb{E} \left(e^{-rT}\max(S_T-K,0)\frac{z}{S\sigma\sqrt{T}} \right)\tag{6}$$

라 쓸 수 있습니다. $z$는 표준정규분포를 따르는 난수이고, $S_T$는 이 난수에서 만들어낸 기초자산 종가이죠.

 

자 그럼 식(6)이 잘 맞는 식인지 코딩을 통해 Black Scholes Formula로 구한 수학 공식과 비교를 해보겠습니다.

 

 

 

Python Code

 

예전 글에서 구현해 놓았던 Closed form formula 관련 코드를 그대로 가져다 쓰겠습니다(관련 코드는 이런 글 등에서 찾아볼 수 있습니다.)

 

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
import pandas as pd

# 아래는 call option delta의 수학 공식
def CallOptionBS(S, K, T, r, q, sigma):
    Ncdf = norm.cdf
    npdf = norm.pdf

    if T == 0:
        val = np.maximum(S - K, 0)
        delta = 1 if S >= K else 0
        gamma = 0
        speed = 0
        theta = 0
        vega = 0
        volga = 0
        ultima = 0

    else:
        d1 = (np.log(S / K) + (r - q + sigma ** 2 / 2) * T) / (sigma * np.sqrt(T))
        d2 = d1 - sigma * np.sqrt(T)
        val = S * np.exp(-q * T) * Ncdf(d1) - K * np.exp(-r * T) * Ncdf(d2)
        delta = np.exp(-q * T) * Ncdf(d1)
        gamma = np.exp(-q * T) * npdf(d1) / (S * sigma * np.sqrt(T))
        speed = -np.exp(-q * T) * npdf(d1) / (S ** 2 * sigma * np.sqrt(T)) * (1 + d1 / (sigma * np.sqrt(T)))
        theta = -np.exp(-q * T) * S * npdf(d1) * sigma / (2 * np.sqrt(T)) - r * K * np.exp(-r * T) * Ncdf(
            d2) + q * S * np.exp(-q * T) * Ncdf(d1)
        # theta2 = -(r - q) * S * delta - 0.5 * sigma ** 2 * S ** 2 * gamma + r * val
        # theta +(r-q)S *delta + 0.5*sigma^2*S^2 *gamma -rf =0
        vega = S * np.exp(-q * T) * npdf(d1) * np.sqrt(T)  # vega
        volga = vega * d1 * d2 / sigma  # volga
        ultima = -vega / sigma ** 2 * (d1 * d2 * (1 - d1 * d2) + d1 ** 2 + d2 ** 2)  # ultima

    # return index
    # 0. value
    # 1. delta
    # 2. gamma
    # 3. speed
    # 4. theta
    # 5. vega
    # 6. volga
    # 7. ultima
    return val, delta, gamma, speed, theta, vega, volga, ultima


def call_option_delta_likelihood(spot, strike, maturity, rfr, div, vol):
    np.random.seed(0)   # 계산떄마다 값이 달라지는 것을 방지하기 위해 seed 고정
    n_iter = 100000     # 10만 개의 sample을 추출해 기댓값을 구할 예정

    drift = rfr - div - 0.5 * vol ** 2   # GBM의 drift
    diffusion = vol * np.sqrt(maturity)  # GBM의 diffusion term
    rand_vec = np.random.normal(size=n_iter)    # standard normal random 100,000 추출
    smat_vec = spot * np.exp(drift + diffusion * rand_vec)  # 만기 종가 생성(10만개)
    
    # 기댓값을 구할 추정량 계산 (식(6)참고)
    estimator = np.maximum(smat_vec - strike, 0) * rand_vec * np.exp(-rfr * maturity) / (spot * vol * np.sqrt(maturity))
    call_delta = estimator.mean()   # 추정량들의 기댓값
    return call_delta


def theoretical_vs_likelihood_call_delta():
    strike = 100
    rfr = 0.02
    div = 0.01
    maturity = 1
    vol = 0.2
    s_vec = np.arange(10, 150, 1)   # 기초자산 가격 변화에 따른 delta 구하려는 의도

    theoretical_delta = [CallOptionBS(s, strike, maturity, rfr, div, vol)[1] for s in s_vec]
    likelihood_delta = [call_option_delta_likelihood(s, strike, maturity, rfr, div, vol) for s in s_vec]
    # 공식 델타와 이 글의 방법론(likelihood)의 delta list
    
    # delta 차 중 가장 많이 차이나는 값 추출
    max_diff = abs(np.array(theoretical_delta) - np.array(likelihood_delta)).max()

    print('maximum error is {:.3f}'.format(max_diff))
    plt.plot(s_vec, theoretical_delta, label='theoretical delta')
    plt.plot(s_vec, likelihood_delta, label='likelihood delta')
    plt.legend()

if __name__ == '__main__':
    theoretical_vs_likelihood_call_delta()

 

 

실행시키면 다음과 같습니다.

 

maximum error is 0.004

 

거의 완전히 일치하는 모양새를 볼 수 있습니다. 델타를 구하는 또 다른  좋은 방법이 되겠네요.

 

아니 근데, 공식이 있으면 공식대로 구하면 되지 이런 복잡한 방법을 알아놓는 게 무슨 장점이 있을까요? 

 

파생상품의 가격 및 델타 등의 공식이 존재하는 경우는 지극히 드뭅니다. 콜옵션은 공식이 존재하는 몇 안 되는 파생상품 중의 하나이죠, 따라서 일반적인 파생상품은 

 

몬테카를로 시뮬레이션

 

방법으로 가격을 구하고,

 

유한차분법

 

으로 델타, 감마 등 민감도를 구합니다(해당 글 참조) 

 

유한차분법으로 민감도를 구하려면 가격을 최소한  2번 구해야 하죠. 일례로 델타 같은 경우에는

$$\frac {c(t, S+\Delta S) - c(t, S-\Delta S)}{2\Delta S}$$

처럼 분자에 위치한 콜옵션 값 2개를 구해 차분을 산출해야 합니다. 

 

반면에 이 글에서 델타 구하는 방법은 가격을 2개 구할 필요 없이 어떤 추정량을 찾아내서 기댓값을 구하는 방법으로도 충분히 민감도를 구할 수 있다는 뜻입니다.

 

따라서 유한차분 방법보다 훨씬 효율적이고 값도 안정성을 띄게 되겠죠.

 

다음 글에서는 감마, 베가 등을 유사한 방법으로 어떻게 구할 수 있는지 알아보겠습니다.