GBM 주가패스 만들기 #1: 만기시점 주가만

 

이 글은

2022.06.07 - [금융공학] - 주식의 수학적 모델 #3 : GBM모델

주식의 수학적 모델 #3 : GBM모델

에서 이어집니다. 주가의 모델을 만들었으니, 실제 주가가 움직이는 흐름, 즉 주가 패스를 어떻게 만들 수 있는지 알아보겠습니다. 주가 패스를 복습하자면, 시점 $t$에서의 주가를 $S_t$라 하고,

• $S_0$ : 현재시점의 주가
• $\mu$ : 주식의 연 수익률
• $\sigma$ : 주식의 연간 수익률의 표준편차, 즉 변동성

이라 할 때,

 

$$ dS_t/S_t = \mu dt + \sigma dW_t $$

가 GBM 모델이라 하였고, 이것을 풀면 미래 시점 $T$에 대해

$$S_T = S_0 \exp\left( \left(\mu-\frac12\sigma^2 \right)T + \sigma W_T \right) \tag{1}$$

또는 미래 시점 수열 $0=t_0 \leq t_1 \leq \cdots \leq t_n=T$ 에 대해

$$ S_{t_{i+1}} = S_{t_i} \exp\left( \left(\mu-\frac12\sigma^2 \right)(t_{i+1}-t_i) + \sigma (W_{t_{i+1}} -W_{t_i}) \right) \tag{2}$$

형태로 recursively 생성할 수도 있다고 하였습니다. 

그럼 python code로 수식(1), (2) 번이 각각 어떤 역할을 하는 식들인지 볼까요?

 

1. 수식(1)처럼 한방에 만기 시점의 $S_T$ 생성하기

두 가지를 볼 것입니다.

1. 첫째는, 만기시점에 예상되는 만기 주가 $S_T$를 한 10,000 개 정도 모아 histogram을 그려 볼 것입니다. 그리고 만기 주가의 평균값이 얼마 나오는지 볼 생각입니다.
2. 둘째로 1번처럼 만기 주가를 10,000개 뽑으면서 만기 주가의 평균이 계속 달라지고 업데이트될 텐데, 그 과정을 보겠습니다.

 

만기시점에 관심이 많은 국군병사들

 

 

python code는 다음과 같습니다.

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

def make_GBM_path():
    s0 = 100
    mu = 0.02
    sigma = 0.2
    T = 1
    nSimulation = 10000
    zNormal = np.random.normal(size=nSimulation)

    drift = (mu - 0.5 * sigma ** 2) * T
    diffusion = sigma * np.sqrt(T)
    s_mat = s0 * np.exp(drift + diffusion * zNormal)
    s_mat_avr = s_mat.mean()
    theoretical_avr = s0*np.exp(mu*T)
    print('Average of S_mat    is {:.4f}'.format(s_mat_avr))
    print('Theoretical average is {:.4f}'.format(theoretical_avr))

    s_step_avr = np.zeros(nSimulation)
    for i in range(nSimulation):
        s_step_avr[i] =s_mat[0:i+1].mean()


    plt.subplot(121)
    plt.hist(s_mat, bins=50)
    plt.axvline(s_mat_avr, color ='r')
    plt.title('historam of $S_T$')

    plt.subplot(122)
    plt.plot(s_step_avr)
    plt.axhline(s_mat_avr, color='r')
    plt.title('Convergence of Average of $S_T$')
    plt.show()

if __name__ == '__main__':
    make_GBM_path()

 

 

간략한 설명은 다음과 같습니다.

    s0 = 100
    mu = 0.02
    sigma = 0.2
    T = 1
    nSimulation = 10000
    zNormal = np.random.normal(size=nSimulation)

현재 시점 주가 $S_0$를 100, 연평균 수익률 $\mu$을 2%, 변동성 $\sigma$를 20%, 만기시점 $T$를 1(year)로 세팅합니다.

이 과정에서 normal random 난수를 nSimulation(=10,000개) 발생시켜 10,000개의 예상 가능한 만기시점 주가를 뽑겠습니다.

 

 

    drift = (mu - 0.5 * sigma ** 2) * T
    diffusion = sigma * np.sqrt(T)
    s_mat = s0 * np.exp(drift + diffusion * zNormal)

수식(1)에서

$$S_T = S_0 \exp\left( \left(\mu-\frac12\sigma^2 \right)T + \sigma \sqrt{T} z \right)~~,~~ z \sim \mathcal{N}(0,1)$$

이고 drift, diffusion term을 각각

$$
\begin{align}
{\rm{drift: }} &= \left(\mu-\frac12\sigma^2 \right)T \\
{\rm{diffustion: }} &=  \sigma \sqrt{T} \\
\end{align}
$$

라 정의하여 식의 계산을 알아보기 쉽게 만든 것입니다.

 

    s_mat_avr = s_mat.mean()
    theoretical_avr = s0*np.exp(mu*T)
    print('Average of S_mat    is {:.4f}'.format(s_mat_avr))
    print('Theoretical average is {:.4f}'.format(theoretical_avr))

10,000개의 만기시점 주가를 얻었다면 그것을 평균하여 s_mat_avr 에 저장합니다.

그런데, 사실 이 평균은 이론적으로 알려져 있거든요. 바로

$$\exp(\mu\cdot T)$$

입니다(Log Normal 분포 블로그를 참고해보세요.) 따라서 이 두 값을 비교해 보는 것입니다.

 

    s_step_avr = np.zeros(nSimulation)
    for i in range(nSimulation):
        s_step_avr[i] =s_mat[0:i+1].mean()

만기 주가를 10,000개 모아놓은 s_mat를 처음부터 $i$개씩 slicing하여 평균을 구합니다. 그러면 만기주가 데이터가 1개에서 10,000개까지 늘어날수록 평균값이 계속 업데이트되겠죠. 상식적으로 처음 몇 개의 평균을 막 흔들리고 이상한 값이 찍히다가 데이터의 개수가 많아질수록 안정적인 어떤 값으로 수렴하게 되겠죠. 그것을 그래프로 보기 위해 s_step_avr라는 변수를 만들고 이것을 ploting 해봅니다.

 

자 그럼 결과를 볼까요?

결과 분석을 해보면 아래와 같습니다.

1. 히스토그램의 경우 만기 주가 $S_T$의 분포는 약간 왼쪽으로 치우친 그럼 모양이 나옵니다. 이것은 수식(1)에서 어느 정도 예상을 할 수 있는데요, 수식(1)의 $S_T$ 공식이 $$\exp(\mathcal{N}(\cdot,\cdot))$$ 형식이라 그런 것입니다. 이런 분포를 log normal 분포라 합니다. 즉, log를 취하면 정규분포이다는 얘기입니다.
2. 오른쪽 그림은 평균이 처음 몇 개의 데이터로는 막 흔들립니다. 80 아래까지 갔다가 110까지 올라갔다가.. 그런데 10,000개까지 데이터가 축척되면서 구해진 평균은 결국 어떤 값으로 안정적으로 수렴하게 되죠. 따라서 어떤 현상의 기댓값을 구할 때는 데이터가 많을수록 안정적인 결과가 나온다라는 것을 말해줍니다.

---

2. 수식(2)처럼 단위 시간으로 쪼개서 주가 패스 생성하기

내용이 길어지므로 다음 글에서 보다 자세하게 다뤄보도록 하겠습니다.