# 選擇權定價-蒙地卡羅模擬法

## 本文重點概要

• 文章難度：★★★★☆
• 使用蒙地卡羅模擬定價選擇權
• 介紹變異數縮減方法以增強定價結果

## 編輯環境與模組需求

``````# Import required packages
import math
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import norm
import time
import tejapi
plt.style.use('bmh')

# Log in TEJ API
api_key = 'YOUR_KEY'
tejapi.ApiConfig.api_key = api_key
tejapi.ApiConfig.ignoretz = True``````

## 資料導入

``````# Import required data
gte, lte = '2021-03-16', '2023-04-20'
stocks = tejapi.get('TWN/APRCD', # stock price
paginate = True,
coid = 'Y9999',
mdate = {'gte':gte, 'lte':lte},
opts = {
'columns':[ 'mdate','close_d']
}
)
# Get options price
puts = tejapi.get( # puts price
'TWN/AOPTION',
paginate = True,
coid = 'TXO202304P15500',
mdate = {'gte':gte, 'lte':lte},
opts = {
'columns':['mdate', 'coid','settle', 'kk', 'theoremp', 'acls', 'ex_price', 'td1y', 'avolt', 'rtime']
}
)
calls = tejapi.get( # calls price
'TWN/AOPTION',
paginate = True,
coid = 'TXO202304C15500',
mdate = {'gte':gte, 'lte':lte},
opts = {
'columns':['mdate', 'coid','settle', 'kk', 'theoremp', 'acls', 'ex_price', 'td1y', 'avolt', 'rtime']
}
)
``````

## 資料處理

``````# Calculate daily return
stocks['daily return'] = np.log(stocks['close_d']) - np.log(stocks['close_d'].shift(1))
stocks['moving volatility'] = stocks['daily return'].rolling(252).std()``````

## 股價預測

1. S0: 初始日之資產價格。
2. r: 資產歷史報酬率。
3. sigma: 資產歷史報酬波動度。
4. T: 到期時間。
5. Nsteps: 切分時間點數量。
6. Nrep: 股價路徑數量。

``````def mc_asset(S0, r, sigma, T, Nsteps, Nrep):
SPATH = np.zeros((Nrep, 1 + Nsteps))
SPATH[:, 0] = S0
dt = T / Nsteps
nudt = (r - 0.5 * sigma **2) * dt
sidt = sigma * np.sqrt(dt)

for i in range(0,Nrep):
for j in range(0,Nsteps):
SPATH[i,j+1] = SPATH[i,j] * np.exp(nudt + sidt * np.random.normal())
return SPATH``````

``````S0 = 100
K = 110
CallOrPut = 'call'
r = 0.03
sigma = 0.25
T = 0.5
Nsteps = 10000
Nrep = 1000
SPATH = mc_asset(S0, r, sigma, T, Nsteps, Nrep)

plt.figure(figsize = (10,8))
for i in range(len(SPATH)):
plt.plot(SPATH[i])
plt.xlabel('Numbers of steps')
plt.ylabel('Stock price')
plt.title('Monte Carlo Simulation for Stock Price')
plt.show()``````

## 選擇權定價

``````def mc_options(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep):
SPATH = mc_asset(S0, r, sigma, T, Nsteps, Nrep)
if CallOrPut == 'call':
payoffs = np.maximum(SPATH[:,-1] - K, 0)
return np.mean(payoffs)*np.exp(-r*T)
else:
payoffs = np.maximum(K - SPATH[:,-1], 0)
return np.mean(payoffs)np.exp(-rT)``````

``````S0 = 100
K = 110
CallOrPut = 'put'
r = 0.03
sigma = 0.25
T = 0.5
Nsteps = 10000
Nrep = 1000
p_ = mc_options(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep)
mybs = BS_formula(S0, K, r, sigma, T)
c, p = mybs.BS_price()

print(f'Monte Carlo price: {c_}')
print(f'Black Scholes price: {p}')``````

## 對偶變數法(Antithetic Variate)

``````def mc_options_AV(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep):

SPATH1 = np.zeros((int(Nrep/2), 1 + Nsteps))
SPATH2 = np.zeros((int(Nrep/2), 1 + Nsteps))
SPATH1[:, 0], SPATH2[:, 0] = S0, S0
dt = T / Nsteps
nudt = (r - 0.5 * sigma **2) * dt
sidt = sigma * np.sqrt(dt)

for i in range(0,int(Nrep/2)):
for j in range(0,Nsteps):
epsilon = np.random.normal()
SPATH1[i,j+1] = SPATH1[i,j] * np.exp(nudt + sidt * epsilon)
SPATH2[i,j+1] = SPATH2[i,j] * np.exp(nudt - sidt * epsilon)

if CallOrPut == 'call':
C1 = np.maximum(SPATH1[:, -1] - K, 0)
C2 = np.maximum(SPATH2[:, -1] - K, 0)
C = np.mean(0.5 * (C1 + C2))
C0 = np.exp(-r*T) * C
return C0
else:
P1 = np.maximum(K - SPATH1[:, -1], 0)
P2 = np.maximum(K - SPATH2[:, -1], 0)
P = np.mean(0.5 * (P1 + P2))
P0 = np.exp(-r*T) * P
return P0``````

``````CallOrPut = 'put'
K = 110
S0 = 100
r = 0.03
sigma = 0.25
T = 0.5
Nrep = 10000
Nsteps = 1000
print('Price under AV: ', mc_options_AV(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep))
print('Price under MC: ', mc_options(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep))``````

## 控制變數法(Control Variate)

• CallOrPut: 選擇權為買或賣權。
• K: 履約價格。
• S0: 期初價格。
• r: 股價歷史報酬。
• sigma: 股價報酬波動度。
• T: 到期時間。
• Nsteps: 切分時間點數量。
• Nrep: 路徑數量。
• Npilot: 計算 c* 所需之時間窗格。
``````def mc_options_CV(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep, Npilot):

# Calculate covariance between stock and options price
SPATH = np.zeros((Npilot, 1 + Nsteps))
SPATH[:, 0] = S0
dt = T / Nsteps
nudt = (r - 0.5 * sigma **2) * dt
sidt = sigma * np.sqrt(dt)

for i in range(0,Npilot):
for j in range(0,Nsteps):
SPATH[i,j+1] = SPATH[i,j] * np.exp(nudt + sidt * np.random.normal())
Sn = SPATH[:, -1]
if CallOrPut == 'call':
Cn = np.maximum(SPATH[:,-1] - K, 0) * np.exp(-r*T)
MatCov = np.cov(Sn, Cn)[0,1]
VarY = S0 ** 2 * np.exp(2 * r * T) * (np.exp(T * sigma ** 2) - 1)
c = -MatCov / VarY
ExpY = S0 * np.exp(r*T)
else:
Pn = np.maximum(K - SPATH[:,-1], 0) * np.exp(-r*T)
MatCov = np.cov(Sn, Pn)[0,1]
VarY = S0 ** 2 * np.exp(2 * r * T) * (np.exp(T * sigma ** 2) - 1)
c = -MatCov / VarY
ExpY = S0 * np.exp(r*T)

# Applying control variate function with optimal c*
SPATH2 = np.zeros((Nrep, 1 + Nsteps))
SPATH2[:, 0] =S0
dt = T / Nsteps
nudt = (r - 0.5 * sigma **2) * dt
sidt = sigma * np.sqrt(dt)

for i in range(0,Nrep):
for j in range(0,Nsteps):
SPATH2[i,j+1] = SPATH2[i,j] * np.exp(nudt + sidt * np.random.normal())
S = SPATH2[:, -1]
if CallOrPut == 'call':
C = np.maximum(SPATH2[:,-1] - K, 0) * np.exp(-r*T)
CVC = np.mean(C + c * (S - ExpY))
return CVC
else:
P = np.maximum(K - SPATH2[:,-1], 0) * np.exp(-r*T)
CVP = np.mean(P + c * (S - ExpY))
return CVP``````

``````CallOrPut = 'put'
K = 110
S0 = 100
r = 0.03
sigma = 0.25
T = 0.5
Nrep = 5000
Nsteps = 1000
Npilot= 5000

print('Price under AV: ', mc_options_AV(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep))
print('Price under CV: ', mc_options_CV(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep, Npilot))``````

## 實際案例

``````S0 = stocks.loc['2023-01-31']['close_d']
K = 15500
r = stocks['daily return'].rolling(252).mean().loc['2023-01-31'] # average return of stock
T = 51 / 252
sigma = stocks.loc['2023-01-31']['moving volatility'] * np.sqrt(252)
Nstep = 50000
Nrep = 50000
Npilot = 5000
CallOrPut = 'put'

print('Monte Carlo theoretical price: ', mc_options(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep))
print('Monte Carlo with AV theoretical price: ', mc_options_AV(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep))
print('Monte Carlo with CV theoretical price: ', mc_options_CV(CallOrPut, K, S0, r, sigma, T, Nsteps, Nrep, Npilot))
print('TEJ Black Scholes price: ', puts.loc['2023-01-31']['theoremp'])
print('Real price: ', puts.loc['2023-01-31']['settle'])``````