GBM모델로 생성된 미래 주가는 어떻게 분포할까? 변동성이 커지면 분포는?

 

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

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

에서 이어집니다. 저번 글에서는 GBM 모델을 따르는 기초자산의 특정 만기 시점 $T$에서의 가격인 $S_T$가 어떤 확률밀도 함수(pdf)를 가지는지, 또 누적분포함수(cdf)는 무엇인지를 알아보았습니다. 잠깐 복습해 보겠습니다.

 

 

복습

기초자산  $S_t$가 GBM 모델 

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

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

 

$S_T$는

$$ \ln(S_T/S_0) = (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_0)-(r-q-\frac12\sigma^2)T}{\sigma\sqrt{T}} $$
확률밀도함수 $g(x)$	$$ g(x)  = \frac{1}{x\sigma\sqrt{T}} \phi(d(x)), $$ 여기서 $\phi(\cdot)$는 표준정규분포의 확률밀도함수

입니다.  누적분포함수와 확률밀도 함수를 그려보면 아래와 같았죠.

 

 

$S_T$의 pdf는 약간 왼쪽으로 쏠린,  소위 log normal 분포로 알려진 그래프를 띄게 됩니다. 

 

그럼 오른쪽 그림의 확률밀도 함수는 변동성의 변화에 따라 어떻게 그려질까요? 상상이 되시는지요.

 

 

 

$S_T$의 변동성과 확률밀도함수간의 관계

 

우선 python coding을 통해 알아보겠습니다.

 

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


def cdf_pdf_of_smat(x, s0, rfr, div, vol, maturity):
    Ncdf = norm.cdf
    npdf = norm.pdf

    drift = (rfr - div - 0.5 * vol ** 2) * maturity
    diffusion = vol * np.sqrt(maturity)

    if x > 0:
        dx = (np.log(x / s0) - drift) / diffusion
        cdf_value = Ncdf(dx)
        pdf_value = 1 / (x * vol * np.sqrt(maturity)) * npdf(dx)
    else:
        cdf_value = 0
        pdf_value = 0

    return cdf_value, pdf_value


def test_of_pdf_vs_volatility():
    s0 = 100     # 기초자산 현재가
    rfr = 0.02   
    div = 0.01
    vol_array = np.arange(0.1, 1 + 0.2, 0.2)   # 변동성 10% ~ 100% 사이의 20% 간격의 변동성
    maturity = 1
    x_vec = np.arange(1, 300, 0.1)   # 1<=기초자산 <=300 으로 변한다고 가정

    max_vec = []              # 각 변동성별로 확률밀도함수의 최댓값이 나오는 x 들 list 
    max_vec_theoretical = []  # 각 변동성별로 확률밀도함수의 이론적 최댓값이 나오는 x list
    for vol in vol_array:
        pdf_vec = []
        for x in x_vec:
            tempvalue = cdf_pdf_of_smat(x, s0, rfr, div, vol, maturity)
            pdf_vec.append(tempvalue[1])  # x 값에 따라, pdf value들 list
        y_max = max(pdf_vec)    # 확률밀도 함수의 최댓값을 구하고
        x_at = x_vec[pdf_vec.index(y_max)]  # 그 최대값이 나오는 index를 찾아 그 index로 x값 찾음

        max_vec.append(x_at)   # 위와 같은 방법으로 찾은 x
        max_vec_theoretical.append(s0 * np.exp((rfr - div - 1.5 * vol ** 2) * maturity))
                               # theretically 찾은 x
        plt.plot(x_vec, pdf_vec, label='vol : {:.0f}%'.format(vol * 100)) # vol vs pdf graph
        plt.scatter(x_at, y_max, s=20)  # 각 pdf의 최댓값에 spotting

    print(max_vec)
    print(max_vec_theoretical)
    plt.legend()
    plt.show()

 

 

우선 그래프를 보면 다음과 같습니다. 변동성이 커질수록,

 

○ 오른쪽 왼쪽으로 더 많이 퍼지는 확률밀도 함수가 됩니다.

 

○ 최댓값을 기준으로 오른쪽 그래프가 완만하게 감소합니다. 그리고 꼬리 부분이 통통 해지죠. 즉, 변동성이 크면 

   주가가 크게 상승할 확률이 높아져서 그렇습니다.

 

○ 최댓값이 나타나는 주가가 더 작아집니다. 변동성이 커질수록 주가가 크게 떨어질 확률이 높지만, 하락의 경우에는 $0$이라는 값이 딱 막고 있죠. 그래서 $S=0$이라는 수직선에 붙어 가는 그래프가 나타나게 됩니다. 

 

또한, 위 코드에서

[99.50000000000009, 88.20000000000007, 69.40000000000006, 48.40000000000004, 30.000000000000025, 16.400000000000013]
[99.50124791926824, 88.24969025845954, 69.41966508779788, 48.43245689553624, 29.969199951324626, 16.447445657715477]

의 출력물을 얻습니다. 

윗 줄은 코딩을 통해 구한 '변동성 별로 pdf의 최댓값이 나타나는 기초자산의 가격을 list 화 시킨 것이고, 

아래 줄은 모종의 수식을 통해 변동성 별로 pdf의 최댓값이 나타나는 기초자산을 차은 list입니다.

 

거의 유사하죠? 값의 차이가 나는 이유는, 코딩에서는

    x_vec = np.arange(1, 300, 0.1)

처럼 discrete 한 $x$들에 대해서 pdf 함수를 구한 것이므로, 이 중에서 최대값을 찾으면 아무래도 정확한 값과 차이를 보이게 되겠죠. 

 

그럼 위 코딩에서 모종의 수식 부분 즉,

        max_vec_theoretical.append(s0 * np.exp((rfr - div - 1.5 * vol ** 2) * maturity))

위 부분의 append로 넣어지는 값인

s0 * np.exp((rfr - div - 1.5 * vol ** 2) * maturity)

이 어떻게 나오는지 아래에서 살펴보겠습니다.

 

 

 

$S_T$ 의 pdf의 최댓값은 어디서 나타날까?

 

위의 표를 보면 $S_T$의 pdf 함수 $g$는 

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

이고, 여기서 함수 $d(x)$는

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

입니다.

 

$g(x)$의 최댓값을 찾기 위해 $g(x)$를 $x$로 미분해 보죠.

$$
\begin{align}
g(x)  &= \frac{-1}{x^2\sigma\sqrt{T}} \phi(d(x)) +\frac{1}{x\sigma\sqrt{T}} \phi'(d(x)) d'(x) \\
      & = \frac{-1}{x^2\sigma\sqrt{T}} \phi(d(x)) -\frac{1}{x\sigma\sqrt{T}} d(x) \phi(d(x)) d'(x) \tag{1}
\end{align}
$$

 

두 번째 식에서 $\phi'(x)=-x \phi(x)$라는 관계식이 쓰였습니다. 또한 

$$d'(x) = \frac{1}{x\sigma\sqrt{T}}$$ 라는 사실과 $g(x)$의 최댓값은

$$ g'(x)=0$$

이라는 관계식에서 구할 수 있으므로 이것들을 식 (1)과 함께 결합해 보면

 

$$ d(x) =-\sigma\sqrt{T}\tag{2}$$

을 만족하는 점 $x$에서 pdf의 최댓값이 나옴을 알 수 있습니다.

 

식(2)을 $x$에 관해 풀면, 

$$ x=S_0 \exp\left((r-q-\frac32\sigma^2)T  \right)$$

를 얻습니다.

 

바로 이것이 code에 나왔던 이론적 공식입니다.