2Star WorstPerform 풋옵션 가격: Closed form #3 (Python Code)

이번 글은

2Star WorstPerform 풋옵션 가격: Closed form #1   

2Star WorstPerform 풋옵션 가격: Closed form #2  

 

에서 살펴보았던 2Star WorstPerform Put option을 직접 코딩해 본 결과를 소개하는 글입니다.

먼저 Closed form을 복기해보죠.

 

 

2Star WorstPerfom Put Option Closed Form

2개의 기초자산 $X_1(t), X_2(t)$ 를 

 

$$ 
\begin{align}
X_i(T) &= X_i(0) \exp\left((r-q_i-\frac12\sigma_i^2)T+\sigma_i W_i(T) \right) \\
       & := x_i \exp (d_i +v_i w_i)~~,~~i=1,2
\end{align}
$$

$$
\begin{align}
x_i & := X_i(0)\\
d_i & := \left(r-q_i =\frac12\sigma_i^2\right)T\\
v_i & := \sigma_i \sqrt{T}
\end{align}
$$

라 세팅했을 때, closed form으로 구한 해당상품의 가격 $p$ 는 다음과 같습니다.

$$p= e^{-rT}[(I)+(II)]$$

$$(I)=  K\Phi_2(\mathbf{p_1},0,D_1) - x_1 \exp \left(d_1+ \frac12 \mathbf{e_1}^t D_1 \mathbf{e_1} \right) \Phi_2(\mathbf{p_1}; D_1\mathbf{e}_1, D_1)  $$
$$(II)= K\Phi_2(\mathbf{p_2},0,D_2) - x_2 \exp \left(d_2+ \frac12 \mathbf{e_2}^t D_2 \mathbf{e_2} \right) \Phi_2(\mathbf{p_2}; D_2\mathbf{e}_2, D_2)  $$ 

기호 정리 ○ $ \mathbf{p_1} =(\ln(K)-\eta_1 , \eta_2-\eta_1), \mathbf{p}_2 = (\eta_1-\eta_2~,~ \ln(K)-\eta_2\} $ ○ $ \Lambda_1  = \begin{pmatrix} v_1 & 0 \\ v_1 & -v_2 \end{pmatrix},  \Lambda_2 = \begin{pmatrix} -v_1 & v_2 \\ 0 &v_2 \end{pmatrix}$ ○ $D_1 :=\Lambda_1 \Sigma \Lambda_1^t , D_2 :=\Lambda_2 \Sigma \Lambda_2^t$  ○ $\mathbf{e}_1 = {\begin{pmatrix} 1 & 0 \end{pmatrix} }^t, \mathbf{e}_2 = {\begin{pmatrix} 0 & 1 \end{pmatrix} }^t$ ○ $d_i = \left(r-q_i-\textstyle{\frac12}\sigma_i^2\right)T~~,~~ v_i = \sigma_i\sqrt{T}~,~i=1,2$ ○ $\Sigma= \begin{pmatrix} 1&\rho\\ \rho &1 \end{pmatrix} $ ○ $\eta_i = \ln(x_i)+d_i~,~i=1,2$

 

이제 코딩해 보겠습니다. 공식이 있으니 그냥 공식에 따라 작성하면 됩니다.

 

 

 

Python Code

 

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

N = norm.cdf

def PutOptionBS(S, K, T, r, q, sigma):
    if T == 0:
        return np.max(K - S, 0)
    else:
        d1 = (np.log(S / K) + (r - q + sigma ** 2 / 2) * T) / (sigma * np.sqrt(T))
        d2 = d1 - sigma * np.sqrt(T)
        return -S * np.exp(-q * T) * N(-d1) + K * np.exp(-r * T) * N(-d2)

# 예전에 다루었던 Market Class, Underlying class
class Market:
    def __init__(self, rfr, correlation):
        self._rfr = rfr
        self._correl = correlation

class Underlying:
    def __init__(self, refprice, spotvalue, volatility, dividend):
        self._volatility = volatility
        self._dividend = dividend
        self._spot = spotvalue
        self._refprice = refprice


from scipy.stats import multivariate_normal as mvn   #다변량 정규분포 계산을 위해

def TwoStarPut_Closedform(objUnderlyings: Underlying, objMarket: Market, strike, maturity):
    nUnderlying = 2
    spot = []
    refprice = []
    vol = []
    div = []
    drift = []
    diffusion = []
    rfr = objMarket._rfr
    corr = objMarket._correl
    corrMatrix = np.array([[1, corr], [corr, 1]])
	
    # 현재가, 기준가, 변동성, 배당률, drift temr(d_i) , diffusion term(v_i)의
    # 벡터화
    for i in range(nUnderlying):
        spot.append(objUnderlyings[i]._spot)
        refprice.append(objUnderlyings[i]._refprice)
        vol.append(objUnderlyings[i]._volatility)
        div.append(objUnderlyings[i]._dividend)
        drift.append((rfr - div[i] - 0.5 * vol[i] ** 2) * maturity)
        diffusion.append((vol[i] * np.sqrt(maturity)))

    spot = np.array(spot)
    refprice = np.array(refprice)
    vol = np.array(vol)
    div = np.array(div)
    drift = np.array(drift)
    diffusion = np.array(diffusion)

    eta = np.log(spot / refprice) + drift                        # 공식의 η_i
    p0 = np.array([np.log(strike) - eta[0], eta[1] - eta[0]])    # 공식의 p1
    p1 = np.array([eta[0] - eta[1], np.log(strike) - eta[1]])    # 공식의 p2
    points = np.array([p0, p1])                                  # (p1, p2) 로 배열화
    lambda0 = np.array([[diffusion[0], 0], [diffusion[0], -diffusion[1]]]) # 공식 Λ1 행렬
    lambda1 = np.array([[-diffusion[0], diffusion[1]], [0, diffusion[1]]]) # 공식 Λ2 행렬

    dMatrix0 = lambda0 @ corrMatrix @ lambda0.transpose()    #공식 중 D1행렬
    dMatrix1 = lambda1 @ corrMatrix @ lambda1.transpose()    #공식 중 D2행렬
    dMatrix = np.array([dMatrix0, dMatrix1])                 # (D1, D1)행렬화
    e0 = np.array([1, 0]).transpose()                        # 공식의 e1
    e1 = np.array([0, 1]).transpose()                        # 공식의 e2
    elementVec = np.array([e0, e1])                          # (e1, e2)배열화

    value = np.zeros(2)      # 공식의 (I)와 (II)를 의미하는 배열 (I, II)

    for i in range(nUnderlying):
        dist1 = mvn(mean=np.zeros(2), cov=dMatrix[i])
        # 평균이 0, covarmatrix가 D_i 인 multivariate normal 분포 설정
        
        dist2 = mvn(mean=dMatrix[i] @ elementVec[i], cov=dMatrix[i])
        # 평균이 D_i e_i , covarmatrix가 D_i 인 multivariate normal 분포 설정
        
        # 위에서 설정한 dist1, dist2 분포의 cdf 함수를 사용하여 계산
        value[i] = strike * dist1.cdf(points[i]) - spot[i] / refprice[i] * np.exp(
            drift[i] + 0.5 * elementVec[i].transpose() @ dMatrix[i] @ elementVec[i]) * dist2.cdf(points[i])

    discountfactor = np.exp(-rfr * maturity)    # 할인팩터

    return discountfactor * (value[0] + value[1])   # e^{-rT}*( (I) + (II) )

 

설명은 주석을 보시면 되겠습니다. 한 가지 설명할 점은 다변량정규분포의 누적분포값을 계산하기 위해  mvn이라는 함수가 쓰였습니다. mvn을 쓰기 위해서는 scipy.stats  library에서 multivariate_normal이라는 함수를 불러와야 합니다. 함수명이 좀 길어서 편하게 쓰기 위해 mvn으로 별명을 짓습니다.

 

from scipy.stats import multivariate_normal as mvn

 

그런 다음

distribution = multivariate_normal(mean = mean, cov= covariance matrix)

처럼 평균벡터  mean과 공분산행렬 covariance matrix를 입력하여 해당 분포를 만들게 됩니다. 이 분포의 이름을 distribution에 정의하고  점 point에서의 확률밀도함숫값(pdf), 누적분포함숫값(cdf)를 계산하기 위해

pdf_value = distribution.pdf(p)
cdf_value = distribution.cdf(p)

이런 식으로 계산하면 됩니다.  scipy.stats.multivariate_normal 문서를 보시기 바랍니다.

 

 

---

 

 

이제 이 함수를 사용하여 두 가지를 관찰해 보도록 하겠습니다.

 

1. 두 기초자산의 상관계수를 -0.9에서 0.9까지 변화시키면서 MC로 구한 값과 위의 closed form code로 구한 값을 비교

 

2. 두 기초자산의 상관계수를 -0.9로 놓고 (나머지 파라미터는 동일) 2star worst put과 1star 풋옵션의 값비교

 

 

if __name__ == '__main__':

    two_assets = []
    underlying_spot1 = Underlying(refprice=100, spotvalue=100, volatility=0.3, dividend=0)
    underlying_spot2 = Underlying(refprice=100, spotvalue=100, volatility=0.3, dividend=0.0)
    # 두 기초자산의 조건을 일치시킴, 기준가=100, 현재가=100, 변동성 30%, 배당률 0%
    # 그리고 상관계수만 다르게 해 볼 예정
    
    rfr = 0.03
    strike = 1
    maturity = 1

    two_assets.append(underlying_spot1)
    two_assets.append(underlying_spot2)

    corr_vec = np.arange(-0.9, 0.9 + 0.1, 0.1)
    put_val_mc_list = []
    put_val_cf_list = []

    for cr in corr_vec:
        market = Market(rfr=rfr, correlation=cr)  #correlation을 -0.9~ 0.9 까지 0.1단위씩
        put_val_mc = TwoStarPut_MC(two_assets, market, strike=strike, maturity=maturity, n_iteration=3 * 10 ** 6)
        put_val_cf = TwoStarPut_Closedform(two_assets, market, strike=strike, maturity=maturity)

        put_val_mc_list.append(put_val_mc)  # correlation 별 mc 값 list
        put_val_cf_list.append(put_val_cf)  # correlation 별 closed form list

    print(np.array(put_val_mc_list) - np.array(put_val_cf_list))  # 두 값의 차이 계산
    # 참고: list끼리 빼기 연산이 안돼, array형으로 변환후 빼기
    
    plt.subplot(1, 2, 1)
    plt.title('MC vs ClosedForm (correlation change')
    plt.plot(corr_vec, put_val_mc_list, label='mc result')
    plt.plot(corr_vec, put_val_cf_list, label='closed form result')
    plt.xlabel('correlation')
    plt.legend()

    put_val_cf_list = []
    put_val_1star_bsf_list = []
    spot_vec = np.arange(50, 150 + 10, 10)   # 이변에는 stock 값을 50 ~ 150 사이에서 변경
    market = Market(rfr=rfr, correlation=-0.999) # 두 기초자산의 상과계수는 -0.999
    # 상관계수을 -1이라 놓지 못하는 이유는 closed form 공식이 오류남

    for spot in spot_vec:
        two_assets = []
        # 이변에도 두 기초자산의 조건을 일치시킴, 기준가=100, 현재가=100, 변동성 30%, 배당률 0%
        underlying_spot1 = Underlying(refprice=100, spotvalue=spot, volatility=0.3, dividend=0)
        underlying_spot2 = Underlying(refprice=100, spotvalue=spot, volatility=0.3, dividend=0.0)
        two_assets.append(underlying_spot1)
        two_assets.append(underlying_spot2)

        put_val_cf = TwoStarPut_Closedform(two_assets, market, strike=strike, maturity=maturity)
        put_val_1star_bsf = PutOptionBS(spot / 100, strike, maturity, rfr, 0, 0.3)

        put_val_cf_list.append(put_val_cf)
        put_val_1star_bsf_list.append(put_val_1star_bsf)

    plt.subplot(1, 2, 2)
    plt.title('1star vs 2star put (correlation =-1)')
    plt.plot(spot_vec, put_val_cf_list, label='2star worst put(corr=-1)')
    plt.plot(spot_vec, put_val_1star_bsf_list, label='1star put BSFormula')
    plt.xlabel('stock')
    plt.legend()
    plt.show()

 

 

결과를 보실까요?

 

 

왼쪽그림 분석

 

○ MC 값과 closed form 값이 완전 겹쳐 보이죠. 어렵게 구한 Closed form이 맞음을 알 수 있습니다. 이렇게 두 가지 이상의 방법으로 다르게 구하여 값이 같은지 보는 것을 크로스 체크라 하고, 반드시 필요한 과정입니다.

 

○ 두 기초자산의 상관계수가 작아질수록 2Star WorstPerform Put 가격이 증가합니다. 그 이유는, 상관계수가 작아질수록 두 기초자산이 반대로 움직인다는 뜻이고, 그렇다면 기초자산의 WorstPerform은 더 하락할 가능성이 높지요, Put 옵션은 하락에 베팅하는 상품으로써, 기초자산 하락 가능성이 높아지기에, 값은 커지는 것입니다.

 

 

오른쪽 그림 분석

 

○  상관계수가 -1이면 두 기초자산이 역조화하여 기초자산이 1개인 put 옵션보다 가치가 커집니다. 그러면 상관계수가 1일 때는 어떨까요? 위의 코드에서 중간의 -0.999 대신 0.999를 대입해 봅니다.

 

market = Market(rfr=rfr, correlation= 0.999)

 

 

와우! 기초자산 1개인 put옵션과 가격이 그냥 일치해 버리죠. 일견 당연한 결과입니다. 상관계수가 1이면 두 기초자산이 똑같이 움직이는 것이므로 그냥 하나의 기초자산이라 볼 수 있기 때문이죠.