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

 

 

 

 

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

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

와 미래주가의 분포를 활용하여 콜옵션 감마 구하기

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

글에 이어 

 

 

콜옵션의 베가

 

를 위 글 들에서 소개했던 방법인 likelihood (우도비) 방법으로 구해보겠습니다. 복습을  다시 한번 해보죠.

 

 

복습

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

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

에서 다룬 내용에서 

기초자산  $S_t$가 GBM 모델 
$$ dS_t/S_t= (r-q) dt +\sigma dW_t$$
($r$는 무위험 이자율, $q$는 연속배당률,  $\sigma$는 변동성(volatility), 그리고 $W_t$는 위너프로세스)

을 따를 때, 만기 시점의 주가 $S_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)$는 표준정규분포의 확률밀도함수

입니다.

 

 

콜옵션의 베가

 

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

$$c(t,S,\sigma)= e^{-rT} \mathbb{E}(\max(S_T-K,0)) \tag{1} $$
입니다.

보통 콜옵션 가격 $c$는 $t, S$에 대한 2변수 함수로 쓰지만, 이 글에서는 변동성 $\sigma$에 대한 민감도 베가(Vega)에 초점을 맞추고 있으므로, $\sigma$까지 포함해서 3 변수 함수로 썼습니다.

 식(1)은

$$
\begin{align}
 c(t,S,\sigma) & = 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}
$$

와 같이 $S_T$의 확률 밀도함수 $g(x)$를 이용하여 쓸 수 있습니다( 변수 $x$는 $S_T$를 의미하고 $x$의 범위가 $(0,\infty)$이므로 적분 구간은 $(0,\infty)$가 됩니다.)

 

콜옵션의 베가는 변동성의 변화에 따른 콜옵션의 변화 즉,

 

$${\rm{Vega}} = \frac{\partial c}{\partial \sigma}$$

이므로 식(2)를 $\sigma$로 편미분해보겠습니다.

 

$$
\begin{align}
\frac{\partial  c}{\partial \sigma} & = \frac{\partial }{\partial \sigma} \int_0^\infty e^{-rT}\max(x-K,0) g(x) dx\\
& =  \int_0^\infty \frac{\partial }{\partial \sigma} \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 \sigma} dx\\
& =  \int_0^\infty e^{-rT}\max(x-K,0) \frac{\partial  g(x)}{\partial \sigma}  \frac{1}{g(x)} \cdot g(x) dx\\
& = \mathbb{E} \left(e^{-rT}\max(x-K,0) \frac{\partial \ln g(x)}{\partial \sigma} \right)\tag{3}
\end{align}
$$

가 됩니다. 이제 위 식 마지막 등식의 수식

$$\frac{\partial \ln g(x)}{\partial \sigma}$$

을 보기좋게 정리해 보겠습니다.

$$g(x) = \frac{1}{x\sigma\sqrt{T}} \phi(d(x))$$
이므로

$$
\begin{align}
\ln g(x) & = \ln \left(\frac{1}{x\sigma\sqrt{T}} \phi(d(x)) \right)\\
& = \ln(\phi(d(x))) -\ln(x\sqrt{T}) -\ln \sigma\\ 
\end{align}
$$

이므로

$$
\begin{align}
\frac{\partial}{\partial \sigma} \ln g(x) &=  \frac{\partial}{\partial \sigma} \ln(\phi(d(x))) -\frac1\sigma\\
& = \frac{\phi'(d(x))}{\phi(d(x))} \frac{\partial d(x)}{\partial \sigma} -\frac1\sigma\\
& = -d(x) \frac{\partial d(x)}{\partial \sigma} -\frac1\sigma~~ (\because \phi'(z)=-z\phi(z))\\
& = -d(x) \frac{\sigma\sqrt{T}-d(x)}{\sigma} - \frac1\sigma \tag{4}
\end{align}
$$

 

마지막 등식은 다음의 Fact에 의한 것입니다.

Fact.

$$ d(x)\sigma\sqrt{T} = \ln(x/S)-(r-q-\frac12\sigma^2)T $$ 의 양변을 $\sigma$로 미분하면 $$\frac{\partial d}{\partial\sigma} \sigma\sqrt{T} + d(x)\sqrt{T} = \sigma {T}$$ 따라서 $$\frac{\partial d}{\partial\sigma} = \frac{\sigma\sqrt{T}-d(x)}{\sigma}$$

입니다.  따라서 식(3)과 식(4)에 의해

 

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

를 만족합니다. 이러한 방법을 likelihood 방법(우도비 방법)이라고 부릅니다.

 

 

 

식(5)은 다음과 같이 이용하면 됩니다.

 

○ 식(8)의 $x$는 만기종가를 뜻합니다. 즉, 어떤 표준정규분포를 따르는 난수 $z$에 대해

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

로 생성할 수 있습니다.

 

○ $d(x)$ 를 계산해 보면,

$$ d(x) = \frac{\ln(x/S)-(r-q-\frac12\sigma^2)T}{\sigma\sqrt{T}} =z$$

 

이죠. 따라서, 

 

 $d(x)$는 표준정규분포 난수, $x$는 이 난수를 통해 만들어지는 GBM 만기 종가

 

라고 이해하면 될 것 같습니다. 이제 코딩을 통해 이 내용을 확인해 보도록 하겠습니다.

 

 

 

 

Python Code

 

예전 글에서 구현해 놓았던 Closed form formula 관련 코드를 그대로 가져다 쓰겠습니다(관련 코드는 이런 글 등에서 찾아볼 수 있습니다.) 또한 아래 코드는 바로 전 글의 코드와 너무나 흡사합니다. 밑의 코드의 estimator 변수만 식(5)의 공식으로 바꿔 주면 됩니다.

 

 

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


# ----------------------------------------
# call option value/Greeks function
# vega, volga, ultima 추가 (theta2 삭제)
# ----------------------------------------
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_vega_likelihood(spot, strike, maturity, rfr, div, vol):
    np.random.seed(0)
    n_iter = 1000000   #simulation 횟수 100만번

    drift = rfr - div - 0.5 * vol ** 2
    diffusion = vol * np.sqrt(maturity)
    rand_vec = np.random.normal(size=n_iter)
    smat_vec = spot * np.exp(drift + diffusion * rand_vec)
    
    # 위의 식(4)의 estimator를 구한다.
    estimator = np.exp(-rfr * maturity) * \
                np.maximum(smat_vec - strike, 0) * \
                (rand_vec ** 2 - rand_vec * vol * np.sqrt(maturity) - 1) / vol
    call_gamma = estimator.mean()
    return call_gamma


def theoretical_vs_likelihood_call_vega():
    strike = 100
    rfr = 0.02
    div = 0.01
    maturity = 1
    vol = 0.2
    s_vec = np.arange(10, 150, 1)

    theoretical_vega = [CallOptionBS(s, strike, maturity, rfr, div, vol)[5] for s in s_vec]
    # closed form 에서 베가는 5th index
    likelihood_vega = [call_option_vega_likelihood(s, strike, maturity, rfr, div, vol) for s in s_vec]

    max_diff = abs(np.array(theoretical_vega) - np.array(likelihood_vega)).max()

    print('maximum error is {:.6f}'.format(max_diff))
    plt.plot(s_vec, theoretical_vega, label='theoretical vega')
    plt.plot(s_vec, likelihood_vega, label='likelihood vega')
    plt.legend()
    plt.show()


if __name__ == '__main__':
    theoretical_vs_likelihood_call_vega()

 

 

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

 

maximum error is 0.104978

감마보다는 덜하지만, 그래도 오차가 작은 편입니다. 그래프는 다음과 같습니다.

 

 

거의 완전히 일치하는 모양새를 볼 수 있습니다. 

 

 

 

다시 한번 말하지만, 유한차분법보다 낫다!

 

만일 MonteCarlo Simulation으로 베가를 구한다고 생각하면

$$ {\rm{Vega}} = \frac{c(t,S,\sigma+\Delta \sigma) - c(t,S,\sigma-\Delta\sigma)}{2\Delta \sigma}$$

와 같이 유한차분법으로 구할 수 있습니다.

 

그렇게 되면 $c(t,S,\sigma+\Delta \sigma)$와 $c(t,S,\sigma-\Delta \sigma)$ 각각을 계산해야 하는 번거로움이 생기고 계산시간도 길어지게 되죠.

 

하지만 위의 우도비 방법을 쓰게 되다면, 훌륭한 추정량(estimator) 하나를 구한 후, 얘의 기댓값을 구하면 되므로 값도 비교적 안정되고 계산 속도도 빨라지는 효과를 누리게 됩니다.