워스트 퍼포머(worst performer)의 분포 #2 : Python Code

이 글에서는

2022.10.26 - [금융공학] - 워스트 퍼포머(worst performer)의 분포 #1

워스트 퍼포머(worst performer)의 분포 #1

에서 이론적으로 다루었던 결과를 몬테카를로 시뮬레이션(Montecarlo Simulation)으로 검증해 보는 내용을 다루어보도록 하겠습니다.

 

MonteCarlo Simulation은 많은 샘플을 뽑아 실험을 행한 결과의 기댓값으로 원하는 값을 산출하는 방법입니다. 구체적인 이론은

2022.08.08 - [금융공학] - Black Scholes Equation의 풀이: 시뮬레이션

Black Scholes Equation의 풀이: 시뮬레이션

에서 확인할 수 있습니다.

 

 

워스트포퍼머 분포 복습

지난 글에서 GBM을 따르는 두 기초자산의 퍼포먼스 $X_1, X_2$가 

$$ dX_i(t)/X_i(t) = (r-q_i)dt +\sigma_i dW_i(t)~,~i=1,2$$

로 모델링 되었을 때,  이 둘의 워스트 포퍼머인 $wp(t) = \min(X_1(t), X_2(t))$의 누적분포함수(cdf) $F_{wp}(K)$는 아래의 식처럼 구할 수 있다고 하였습니다($x_1, x_2$는 두 퍼포먼스의 초기 가격입니다.)

기초자산 $X_1, X_2$ 의 워스트포퍼머는 다음과 같은 누적분포함수 $F_{wp}(K)$를 따른다.

$$F_{wp}(K) = \Phi(\alpha_1)+\Phi(\alpha_2) -\Phi_2(\alpha_1, \alpha_2 ; \rho)$$
여기서
$$ \alpha_i = \frac{\ln(K/x_i)-\mu_i t}{\sigma_i\sqrt{t}}~,~i=1,2 \tag{1}$$
이다.

 

그럼 이것을 Montecarlo로 어떻게 구할 수 있을까요?

 

 

시뮬레이션 절차

1. 두 기초자산 퍼포먼스 생성하기

기초자산 퍼포먼스 $X_1 , X_2$는  

$$ X_i(t) = x_i \exp \left( (r-q_i -\textstyle{\frac12}\sigma_i^2) t +\sigma_i \sqrt{t} z_i\right)~~,i=1,2$$

의 형태로 생성할 수 있습니다. 여기서 $z_i \sim \mathcal{N}(0,1)$이고 $z_1$과 $z_2$의 상관계수는 $\rho$입니다.

 

$z_1$과 $z_2$에 상관관계를 주기 위해,  [금융공학] - 상관관계를 보이는 두 자산의 움직임 모델 에서 살펴본 대로

○ 서로 독립인 표준 정규분포 난수 $w_1, w_2$를 추출하고

○ $z_1 = w_1~,~ z_2=\rho w_1 + \sqrt{1-\rho^2}w_2$ (촐레스키 분해)

로 $z_1, z_2$를 생성하면 됩니다.

 

이제, 적절한 시뮬레이션 횟수(예를 들어 $N=10,000$회)를 설정하여 $(z_1, z_2)$ 페어를 시뮬레이션 횟수 개만큼 뽑습니다.

 

 

2. 워스트 포퍼머 및 분포 구하기

1번 과정에서 뽑은 $N$개의 퍼포먼스 페어

$$ (X_1^{(1)}, X_2^{(1)}) , \cdots , (X_1^{(N)}, X_2^{(N)}) $$

에 대해 각각의 워스트 포퍼머를 추출해 냅니다.

$$ wp^{(i)} :=\min\left(X_1^{(i)}, X_2^{(i)}\right)~~,~~i=1,\cdots, N $$

 

이제 고정된 값 $K$에 대해 

$$ wp^{(i)} <K$$ 인 것의 개수를 구해서 전체 횟수 $N$ 회 중의 비율을 구하면 바로 $K$에서의 누적분포함숫값이 됩니다.

이것을 위의 식(1)과 비교해 보면 됩니다.

 

 

 

Python Code

import numpy as np
from scipy.stats import multivariate_normal
from scipy.stats import norm


def WorstPerformer_test():
    spot = np.array([1, 1])
    q = np.array([0, 0])
    vol = np.array([0.2, 0.3])
    rfr = 0.02
    mat = 1
    rho = 0.5

    strike = 0.8

    drift = (rfr - q - 0.5 * vol * vol) * mat
    diffusion = vol * np.sqrt(mat)

    norm_rand = np.random.normal(size=(10000, 2))
    norm_rand_correlated = norm_rand
    norm_rand_correlated[:, 1] = rho * norm_rand[:, 0] + np.sqrt(1 - rho ** 2) * norm_rand[:, 1]

    mat_price = spot * np.exp(drift + diffusion * norm_rand_correlated)
    res = np.min(mat_price, axis=1)
    res_below = res[res < strike]
    print('Montecarlo Simulation Result')
    print(len(res_below) / len(res))

    alpha = (np.log(strike / spot) - (rfr - q - 0.5 * vol * vol) * mat) / vol / np.sqrt(mat)

    mul_mean = np.zeros(2)
    mul_cov = np.identity(2)
    mul_cov[0, 1] = mul_cov[1, 0] = rho

    multinorm = multivariate_normal(mul_mean, mul_cov)
    closedform_value = norm.cdf(alpha[0]) + norm.cdf(alpha[1]) - multinorm.cdf(alpha)
    print('Cloase Form Result')
    print(closedform_value)


if __name__ == '__main__':
    WorstPerformer_test()

 

 간략히 살펴보겠습니다.

 

    spot = np.array([1, 1])	# 초기값
    q = np.array([0, 0])	# 배당률
    vol = np.array([0.2, 0.3])	# 변동성
    rfr = 0.02	#무위험 이자율
    mat = 1	# 잔존만기
    rho = 0.5	# 상관계수
	# 퍼포먼스 모델링을 위한 parameter 설정

 

---

    drift = (rfr - q - 0.5 * vol * vol) * mat	# 기초자산 #1,#2의 drift term
    diffusion = vol * np.sqrt(mat)	# 기초자산 #1,2의 diffusion term

 

---

    norm_rand = np.random.normal(size=(10000, 2))	# N=10,000개의 서로 독립인 표준정규난수 (z1,z2) 추출
    norm_rand_correlated = norm_rand
    norm_rand_correlated[:, 1] = rho * norm_rand[:, 0] + np.sqrt(1 - rho ** 2) * norm_rand[:, 1]
    # Choelsy 분해를 통해 서로독립인 (z1,z2)를 상관관계가 있는 페어 (w1,w2)로 변경

 

---

    mat_price = spot * np.exp(drift + diffusion * norm_rand_correlated)	
    # 기초자산1,2 종가생성 (X1,X2), mat_price는 10,000*2 사이지의 배열
    res = np.min(mat_price, axis=1)	# 열방향으로, 즉, (X1,X2) 의 최솟값을 구함, 10,000*1 배열
    res_below = res[res < strike]	# 값 strike보다 워스트포퍼머가 밑에 있는 배열 추출
    print('Montecarlo Simulation Result')
    print(len(res_below) / len(res))	# (wp<strike)개수를 N으로 나눠서 확률 구함

 

---

    alpha = (np.log(strike / spot) - (rfr - q - 0.5 * vol * vol) * mat) / vol / np.sqrt(mat)
		# 식(1)의 alpha1. alpha2 구함
        
    mul_mean = np.zeros(2)	# 표준정규분포를 위해 평균은 0벡터, 
    mul_cov = np.identity(2)	# 각각의 표준편차도 1로 하고, 상관계수만 살림
    mul_cov[0, 1] = mul_cov[1, 0] = rho

○ 평균은 0 벡터, 공분산 행렬은 $\begin{pmatrix} 1&\rho\\ \rho&1 \end{pmatrix}$ 로 합니다. 

 

---

    multinorm = multivariate_normal(mul_mean, mul_cov)
    closedform_value = norm.cdf(alpha[0]) + norm.cdf(alpha[1]) - multinorm.cdf(alpha)
    # 식(1) 사용
    print('Cloase Form Result')
    print(closedform_value)

○ scipy.stats에 있는 multivariate_normal 클래스는 다변량 정규분포에 관련된 정보를 제공하는 함수입니다. 파라미터는 평균 벡터와 공분산을 받습니다.

○ 다변량 정규분포의 누적분포함수를 구하기 위해서는 multivariate_normal.cdf() 를 사용하면 됩니다.

 

 

 

결과를 보실까요?

 

Montecarlo Simulation Result
0.3168
Cloase Form Result
0.3120016862157784

Process finished with exit code 0

seed를 고정하지 않아 실행할 때마다 MC 값은 다르게 나오지만, 식(1)의 결과와 거의 같게 나옴을 알 수 있습니다.

 

 

 

이제 식(1)을 사용할 수 있다는 안도감이 좀 생기네요. 이렇게 여러 가능한 방법으로 교차 검증이 되면 값에 확신이 생기게 됩니다.