스터디/Tech

[기술서적] 파이썬을 활용한 금융 분석 - 확률 과정

_leezoee_ 2024. 3. 10. 14:07

 

 

 

난수 생성

시뮬레이션 작업의 근간이 되는 난수 생성 작업.

Sobol 수열에 기반한 준-난수, 일반적인 의사 난수(Pseudo-random number) 등이 있다.

 

 

import math
import numpy as np
import numpy.random as npr
from pylab import plt, mpl

plt.style.use('seaborn')
mpl.rcParams['font.family'] = 'serif'
%matplotlib inline

# 시드값 고정 및 출력 자릿수 설정
npr.seed(100) 
np.set_printoptions(precision=4)

# 균일분포 난수 1차원 ndarray 객체
npr.rand(10)
# 균일분포 난수 2차원 ndarray 객체
npr.rand(5, 5)

a = 5. #하한
b = 10. #상한
npr.rand(10) * (b - a) + a  #구간 변환
npr.rand(5, 5) * (b - a) + a #2차원 변환

 

 

표준정규분포는 분석이나 수치해석 분야에서 필수불가결한 도구고 가장 많이 사용된다.

네가지 예시 구현.

 

1. 평균이 0이고 표준편차가 1인 표준정규분포

2. 평균이 100이고 표준편차가 20인 정규분포

3. 0.5 자유도를 가진 카이제곱분포

4. 람다계수가 1인 포아송 분포

 

 

히스토그램 시각화 진행

sample_size = 500
rn1 = npr.standard_normal(sample_size)
rn2 = npr.normal(100, 20, sample_size)
rn3 = npr.chisquare(df=0.5, size=sample_size)
rn4 = npr.poisson(lam=1.0, size=sample_size)

fig,((ax1, ax2) , (ax3, ax4)) = plt.subplots(nrows=2, ncols=2, figsize=(10, 8))

ax1.hist(rn1, bins=25)
ax1.set_title('standard normal')
ax1.set_ylabel('frequency')
ax2.hist(rn2, bins=25)
ax2.set_title('normal(100, 20)')
ax3.hist(rn3, bins=25)
ax3.set_title('chi square')
ax3.set_ylabel('frequency')
ax4.hist(rn4, bins=25)
ax4.set_title('Poisson')

 

네가지 분포에 대한 난수 히스토그램

 

 

 

 

 

시뮬레이션

몬테카를로 시뮬레이션은 금융공학에서 가장 중요하고 널리 쓰이는 수치해석 방법 중 하나.

몬테카를로 시뮬레이션이 적분 등의 수학식이나 금융 파생상품의 가치를 계산할 때 가장 유연한 수치해석 방법론이기 때문.

 

 

* 확률 변수

- 옵션 가격 평가를 위한 '블랙-숄즈-머튼 모형'

블랙-숄즈-머튼 모형을 사용하여 미래 특정 시점 T에서의 주가 를 계산하는 공식

 

시뮬레이션 코드

S0 = 100 # 현재 주가
r = 0.05 # 상수 무위험 단기 이자율
sigma = 0.25 # 상수 변동성
T = 2.0 # 시간 (연 단위 만기)
I = 10000 # 시뮬레이션 횟수

#미래시점 T에서의 주가 저장, 이산화 방식은 npr.standard_normal() 사용
ST1 = S0 * np.exp((r - 0.5 * sigma ** 2) * T + 
                 sigma * math.sqrt(T) * npr.standard_normal(I))

plt.figure(figsize = (10,6))
plt.hist(ST1, bins=50)
plt.xlabel('index level')
plt.ylabel('frequency')

 

standard_normal 블랙-숄즈-머튼 모형 주가 시뮬레이션, 기하 브라운 운동 모형

 

 

#npr.lognormal() 함수 사용
ST2 = S0 * npr.lognormal((r - 0.5 * sigma ** 2) * T,
                        sigma * math.sqrt(T) , size=I)

plt.figure(figsize = (10,6))
plt.hist(ST2, bins=50)
plt.xlabel('index level')
plt.ylabel('fequency')

 

standard_lognormal 블랙-숄즈-머튼 모형 주가 시뮬레이션, 기하 브라운 운동 모형

 

import scipy.stats as scs

def print_statistics(a1, a2):
    sta1 = scs.describe(a1)
    sta2 = scs.describe(a2)
    print('%14s %14s %14s' %
         ('statistic' , 'data set 1' , 'data set 2'))
    print(45 * "-")
    print('%14s %14.3f %14.3f' % ('size' , sta1[0], sta2[0]))
    print('%14s %14.3f %14.3f' % ('min', sta1[1][0], sta2[1][0]))
    print('%14s %14.3f %14.3f' % ('max', sta1[1][1], sta2[1][1]))
    print('%14s %14.3f %14.3f' % ('mean', sta1[2], sta2[2]))
    print('%14s %14.3f %14.3f' % ('std', np.sqrt(sta1[3]), 
                                  np.sqrt(sta2[3])))
    print('%14s %14.3f %14.3f' % ('skew', sta1[4], sta2[4]))
    print('%14s %14.3f %14.3f' % ('kurtosis', sta1[5], sta2[5]))
    
print_statistics(ST1, ST2)

두 시뮬레이션 결과 통계치 비교

 

 

 

* 확률 과정

 

확률 변수의 수열, 확률 과정을 시뮬레이션 한다는 것은 확률 변수를 반복해서 시뮬레이션하는 것과 유사.

각 단계의 샘플이 독립적이지 않고 이전 결과들에 의존한다는 차이점이 있음.

금융에서 사용되는 대부분 확률 과정은 미래 값이 바로 전 과거 값에만 의존하고 그보다 더 이전이나 전체 경로에는 의존하지 않는 '마코프 특성'을 가지는 확률 과정이다. 이러한 확률 과정을 '무기억성 확률 과정' 이라고 함.

 

- 기하 브라운 운동 모형

 

블랙-숄즈-머튼 모형에서 사용되는 확률 과정 중 하나로서, 로그-정규 분포를 따름

 

위 수식의 확률 미분 방정식은 오일러 방식으로 정확히 이산화 할 수 있다.

 

I = 10000 #시뮬레이션 반복 횟수
M = 50 #시간 구간을 나누는 수
dt = T/M #시간 구간의 길이, 전체 시간 T를 50으로 나눈 값
S = np.zeros((M + 1, I)) #주가를 저장하는 배열 생성, 각 열은 주가 시뮬레이션 경로를 나타냄
S[0] = S0 #초기 주가를 설정

#현재 시간에서의 주가를 이전 시간에서의 주가에 주가 변화량을 더해 계산, 주가 변화량은 블랙-숄즈-머튼 모형을 따름.
for t in range(1, M+1) : #시간 구간을 순회
    S[t] = S[t-1] * np.exp((r - 0.5 * sigma ** 2) * dt + 
                          sigma * math.sqrt(dt) * npr.standard_normal(I))

plt.figure(figsize = (10,6))
plt.hist(S[-1], bins=50)
plt.xlabel('index level')
plt.ylabel('frequency')

 

블랙-숄즈-머튼 모형에 따라 동적인 주가 시뮬레이션 수행

 

print_statistics(S[-1], ST2)

동적, 정적 시뮬레이션 결과 통계치 비교

 

#블랙-숄즈-머튼 모형에 따라 시뮬레이션 된 주가 경로를 시간에 따라 그래프로 시각화
plt.figure(figsize = (10,6))
plt.plot(S[: , :10], lw=1.5)
plt.xlabel('time')
plt.ylabel('index level')

동적 시뮬레이션 주가 경로 시각화

 

 

- 제곱근 확산 모형

 

금융공학에서 단기 이자율, 변동성 모형에 사용되는 '평균 회귀 과정'

평균 회귀 과정에서 가장 널리 사용되는 모형은 콕스-잉거솔-로스의 제곱근 확산 모형.

이자율 확산 모형

 

이 제곱근 확산 모형을 계산상 장점이 많은 오일러 방식을 사용해 이산화한다.

이산화를 통해 시간에 따른 주가 변화를 계산하고 모의할 수 있다.

 

제곱근 확산 모형은 Xt 값이 항상 양수로 유지되는 편리한 특성이 있으나, 단순 오일러 방식으로 이산화 시 음수가 될 수 있음. 이를 방지하고자 원래 시뮬레이션값에서 양수만 취하는 방식을 사용.

시뮬레이션 코드에 확률과정의 값을 저장하기 위해 하나가 아닌 두 개의 ndarray 객체 사용.

 

x0 = 0.05 #단기 이자율 초깃값, 첫번째 시간에서의 이자율
kappa = 3.0 #평균 회귀계수(회귀 속도)
theta = 0.02 #장기 평균 회귀 값
sigma = 0.1 #변동성
I = 10000 #시뮬레이션 반복횟수
M = 50 #시간 구간의 수
dt = T/M #각 시간 구간의 길이

def srd_euler():
    xh = np.zeros((M+1, I)) #이산화된 이자율 값을 저장하는 배열
    xh[0] = x0
    for t in range(1, M+1):
        xh[t] = (xh[t-1] +
                 kappa * (theta - np.maximum(xh[t-1], 0)) * dt +
                 sigma * np.sqrt(np.maximum(xh[t-1], 0)) *
                 math.sqrt(dt) * npr.standard_normal(I))
    x = np.maximum(xh, 0)
    return x

x1 = srd_euler()

plt.figure(figsize=(10, 6))
plt.hist(x1[-1], bins=50)
plt.xlabel('value')
plt.ylabel('frequency')

오일러 방식을 적용한 제곱근 확산 모형 시뮬레이션 결과 만기 값 시각화

 

 

이 확률 과정을 완벽하게 이산화 시뮬레이션 하기 위해서는

자유도와 비중심 인수를 가진 비중심 카이 제곱 분포로 표현해야한다. 

"""
정학한 해를 이용해 제곱근 확산 모형에 따라 주어진 매개 변수들을 사용해 시뮬레이션 수행.
주어진 시간 간격 내에서 다음 과정을 거쳐 이자율 값을 갱신
1. 자유도와 비중심인수를 계산
2. 비중심 카이제곱분포 난수 생성해 주가 갱신
계산된 주가는 x 배열에 저장
"""

def srd_exact():
    x = np.zeros((M + 1, I))
    x[0] = x0
    for t in range(1, M+1) :
        df = 4 * theta * kappa / sigma ** 2
        c = (sigma ** 2 * (1 - np.exp(-kappa * dt))) / (4 * kappa)
        nc = np.exp(-kappa * dt) / c * x[t-1]
        x[t] = c* npr.noncentral_chisquare(df, nc, size=I)
    return x
    
x2 = srd_exact()   

plt.figure(figsize = (10,6))
plt.hist(x2[-1] , bins=50)
plt.xlabel('value')
plt.ylabel('frequency')

 

정확한 이산 방식을 적용한 제곱근 확산 모형 시뮬레이션 결과 만기 값 시각화

print_statistics(x1[-1] , x2[-1])

오일러 이산화, 정확한 이산화 방식 성능 비교

 

두가지 방식의 통계치를 비교하면 큰 차이는 없음.

I = 250000
%time x1 = srd_euler()
%time x2 = srd_exact()

두 방식 계산 시간 비교

 

계산 시간은 정확한 이산화 방법이 오일러 방식보다 약 두배가 소요된다.

 

 

- 확률적 변동성 모형

 

주식의 변동성이 시간이 지남에 따라 변화할 수 있다는 가정을 반영.

상관계수를 통해 시장이 하락할 때 변동성이 증가하고 시장이 상승할 때 변동성이 감소하는 '레버리지 효과'를 설명.

 

1. 상관계수 행렬 : 각 확률 과정의 상호간 상관계수를 행렬 형태로 나타냄

2. 숄레스키 분해 : 상관계수 행렬을 숄레스키 분해하여 얻은 하삼각행렬은 각 확률 과정 간 상관관계를 나타냄. 이 하삼각행렬은 각 행렬 요소가 각 확률 과정의 상관계수를 나타내며, 이를 통해 두 확률 과정 간 상관관계를 직접적으로 알 수 있음.

#주어진 상관계수 행렬을 이용해 숄레스키 분해를 수행

S0 = 100. #초기 주가
r= 0.05 #이자율
v0 = 0.1 #초기 변동성
kappa = 3.0 #회귀속도
theta = 0.25 #장기평균 회귀값
sigma = 0.1 #변동성
rho = 0.6 #두 확률 과정 간 상관계수 
T = 1.0 #만기 시간

corr_mat = np.zeros((2,2))
corr_mat[0, :] = [1.0, rho]
corr_mat[1, :] = [rho, 1.0]
cho_mat = np.linalg.cholesky(corr_mat) #행렬의 숄레스키 분해

cho_mat
#제곱근 확산 모형을 사용해 변동성을 시뮬레이션 하고, 시간에 따른 변동성의 추정치를 계산하는 과정

M = 50
I = 10000
dt = T/M

ran_num = npr.standard_normal((2, M+1, I)) #3차원 난수 데이터 생성

v = np.zeros_like(ran_num[0])
vh = np.zeros_like(v)

v[0] = v0
vh[0] = v0

for t in range(1, M+1):
    ran = np.dot(cho_mat, ran_num[: , t, :]) #고나련 난수만 고르고 숄레스키 행렬로 변환
    vh[t] = (vh[t-1] +
            kappa * (theta - np.maximum(vh[t-1] , 0)) * dt +
            sigma * np.sqrt(np.maximum(vh[t-1] , 0)) * 
            math.sqrt(dt) * ran[1]) #오일러 방식을 사용한 경로 시뮬레이션
    
v = np.maximum(vh, 0)
#주가 과정도 상관관계를 고려하고 오일러 방식으로 이산화 진행

S = np.zeros_like(ran_num[0]) #주가 배열 초기화
S[0] = S0 #초기 주가 값

#변동성과 주가를 시간에 따라 업데이트, 시간 단계마다 변동성에 대한 확률과정을 시뮬레이션하고 이를 사용해 주가 업데이트
for t in range(1, M+1):
    ran = np.dot(cho_mat, ran_num[: , t, :])
    S[t] = S[t-1] * np.exp((r - 0.5 * v[t]) * dt +
                           np.sqrt(v[t]) * ran[0] * np.sqrt(dt))
    
#시각화
fig, (ax1, ax2) = plt.subplots(1, 2, figsize = (10,6))
ax1.hist(S[-1] , bins = 50)
ax1.set_xlabel('index label')
ax1.set_ylabel('frequency')
ax2.hist(v[-1], bins=50)
ax2.set_xlabel('volatility')

주가, 변동성 시뮬레이션 결과 시각

 

오일러 방식을 사용하면 표준정규분포를 사용하므로 샘플 생성 시 상관관계를 간단하고 일괄적으로 적용할 수 있다.

예시로 주가 과정에는 오일러 방식, 변동성 과정에는 비중심 카이제곱 방식 이러면 상관관계를 일관성 있게 적용하기 어렵다.

 

 

- 점프확산

 

시장에서 볼 수 있는 특징으로 자산 가격이나 변동성이 점프하는 점프확산 현상이 있다.

로그 정규분포에 점프 요인을 더해 블랙-숄즈-머튼 모형을 개선한 점프확한 모형을 발표.

 

점프 확산 모형을 시뮬레이션하기 위해 독립적인 세 개의 난수 집합이 필요

점프로 인한 두 개 최고점 분포를 확인할 수 있는 시뮬레이션 코드를 진행하였다.

 

S0 = 100.  
r = 0.05
sigma = 0.2
lamb = 0.75 #점프강도
mu = -0.6 #평균 점프 크기
delta = 0.25 #점프 변동성
rj = lamb * (math.exp(mu + 0.5 * delta ** 2) - 1) #표류계수 수정치

T = 1.0
M = 50
I = 10000
dt = T / M

S = np.zeros((M + 1, I))
S[0] = S0
sn1 = npr.standard_normal((M  + 1, I)) #표준정규분포 난수
sn2 = npr.standard_normal((M + 1 , I)) #표준정규분포 난수
poi = npr.poisson(lamb * dt, (M + 1, I)) #포아송 분포 난수

for t in range(1, M + 1, 1): #정확한 오일러 방식에 따른 시뮬레이션
    S[t] = S[t-1] * (np.exp((r - rj - 0.5 * sigma ** 2) * dt +
                           sigma * math.sqrt(dt) * sn1[t]) + 
                    (np.exp(mu + delta * sn2[t]) - 1) *
                    poi[t])
    S[t] = np.maximum(S[t] , 0)
    

plt.figure(figsize = (10, 6))
plt.hist(S[-1] , bins = 50)
plt.xlabel('value')
plt.ylabel('frequency')

 

점프확산 모형 시뮬레이션 결과 시각화