Table of Contents
蒙地卡羅模擬法被廣泛運用於財務金融領域，股價路徑預測 中我們已經介紹如何使用蒙地卡羅模擬預測股票價格，今日我們將推廣到更難的選擇權定價模型。在 CRR模型 與 Black-Scholes 模型與 Greeks 兩篇文章中，我們分別使用二元樹原理與 Black-Scholes 公式解計算買賣權理論價格，其中也解釋了許多選擇權相關的基礎知識，若對選擇權一知半解的小夥伴，建議先閱讀完這兩篇文章後，再回來閱讀本篇文章。於之後的章節，我們會先演示如何使用蒙地卡羅法預測股價，再來推廣到歐式選擇權定價，最後介紹一些變異數縮減方法輔助我們使用蒙地卡羅模擬。
本文使用 Window 作業系統以及 Jupyter Notebook 作為編輯器。
# 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公司交易面資料庫: 未調整股價(日)，資料代碼為(TWN/APRCD)。
衍生性金融商品資料庫: 選擇權日交易狀況，資料代碼為(TWN/AOPTION)。
使用台灣加權股價指數(Y9999)未調整收盤價，時間區間為2021/01/31到2023/04/19。並且載入台灣加權指數買權與賣權(TXO202304C15500、TXO202304P15500)，該選擇權為歐式買賣權、開始交易日為1/31，到期日為4/19，履約價格為15500，並且將 “mdate” (日期)欄位設定成索引。
# 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']
}
)
計算大盤之日報酬並且計算移動報酬波動度，以252天也就是一年為窗格迭代下去。
# 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()我們使用蒙地卡羅模擬股價路徑。蒙地卡羅模擬法概念十分簡單，就是先取得資產的報酬過程並且將其離散化後，以極小時間區段分區推算資產價格變化。例如以股票價格為例，其報酬服從幾何布朗運動，因此可以得到其離散化隨機微分過程(式一)，其中的 Wt 為維納過程。再經過伊藤公式後，就可以得到式二作為蒙地卡羅模擬法的主要方程式進行股價預測，其中的 Zt 為標準常態分配。
再來我們可以使用 Python 將上式程式化，蒙地卡羅法精髓在於我們同時以式 2 估計多個股價路徑，最後將每個路徑的最後一個股價加總平均，就會得到是我們欲預測之股價。這裡以我們定義以下變數:
並且寫成 “mc_asset” 函式以執行蒙地卡羅模擬，程式碼如下所式。
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完成函式編輯後，我們帶入數字檢驗，並且繪製結果成圖 1，圖 1 每條線就是一個股價過程。
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)其中我們可以使用 CallOrPut 參數調整選擇權種類，我們帶入數值後進行驗算，並且比較使用 Black-Scholes 模型與 Greeks 所使用之方法下，所得理論價格是否有出路。
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}')所得結果如下圖 2 所示，可以發現兩者確實有些許差距但不大。
由於蒙地卡羅模擬選擇權價格就類似於生成大量可能的理論價格再進行平均，這樣做必定會受到模擬出來價格波動度較大的問題，也就是模擬出來的價格可能會出現許多極端值，進而造成模擬結果失準。為了解決這項問題，我們可以採用對偶變數法降低波動度。其觀念就是當我們生成一條股價路徑(式 3)時，同時生成一個反向的報酬路徑(式 4)，這時候兩條路徑的相關係數為 -1，造成兩者相加所得的共變異數為最小，進而降低預估出選擇權價格的波動度。
具體程式碼編寫如下，我們先生成兩個矩陣進行計算，SPATH1 為第一個矩陣且為正向，SPATH2 則是反向所組成之矩陣，可以發現每次迭代生成隨機亂數 (epsilon) 時，兩個矩陣皆共享亂數，因此所需計算量降低並且降低選擇權預測價格之波動程度。
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接著代入數值檢驗在對偶變數法之下的選擇權價格與一般蒙地卡羅之結果是否一致，結果請見圖 3，可以得出確實十分接近。
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))
除了對偶變數法以外，也可以透過控制變數法來降低選擇權理論價格的波動程度。假設有一兩隨機變數 X 與 Y，其中 Y 變數平均數與變異數計算上十分容易。又假設兩變數可以組成新變數 Z (式 5)，此時 Z 的期望值與 X 期望值相同(見式 6)，而變異數則由 c 決定，因此可以找出一個最佳的 c* 使 Z 的變異數達到最小，c* 最佳解如式 7。我們將 X、Y 視為每個路徑上選擇權與股票現貨價格，我們使用過去選擇權與股價的共變異數與股價之變異數計算最佳的 c*，接著再透過式 5 計算選擇權理論價格 (E(Z))，就可以達成縮減蒙地卡羅波動度的目的。
程式化結果如下方所表示:
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透過代入數值，我們檢驗控制變數法與對偶變數法之下，所得的賣權價格結果如圖 4，可發現幾乎相同。
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))
總結上述，我們學到了總共三種使用蒙地卡羅模擬計算選擇權理論價格的方法，分別為: 一般法、對偶變數法與控制變數法，接著我們就導入真實案例，對比各方法的理論價格與 TEJ 所計算的 Black Sholes 價格是否有很大的差異。
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'])我們採用時間為 2023-01-1 到 2023-04-19、履約價格為 15500 的台指賣權，Sigma 採用以 252 天為窗格的台指報酬率標準差，r 為最近 252 天的平均報酬率，且假設今日 2023 年 1 月 31 日。所得結果如下圖 5，可以發現三種方法所得結果與 TEJ 計算的價格相比，都更趨近於實際價格。
蒙地卡羅定價法相較於 CRR 模型與 Black-Scholes 模型，更仰賴大數法則，通過大量的股價路徑，逐漸擬合到合理的理論價格。在這個電腦效能大幅進步的今天，這種資料驅動的演算法、定價模型相信只會越來越多。投資人在進行選擇權投資時，不彷也將蒙地卡羅法也納入考量。之後我們會提供更多選擇權或衍生性商品相關的知識，歡迎持續關注本平台，此外，也歡迎讀者及投資者們選購 TEJ E Shop 中的方案來建構自己的選擇權定價程式。請注意以上訂價公式與選擇權商品僅為示範所用，不代表本台任何投資或標的上的建議。
