兩種資產最大值的看漲期權的蒙地卡羅近似
我想用收益計算期權的價格 $$ \begin{equation} \max \big{\max{S^1_T, S^2_T} - K, 0\big}, \end{equation} $$ 在哪裡 $ S^{1,2} $ 具有相同的動態，相關性為 0。所以， $$ \begin{align} dS^1_t &= r S_t^1 dt + \sigma S^1_t dW^1_t \ dS^2_t &= r S_t^2 dt + \sigma S^2_t dW^2_t, \end{align} $$ 在哪裡 $ W^1 $ 和 $ W^2 $ 是定價度量下的獨立維納過程 $ Q $ . 該期權具有分析定價公式（例如，在期權定價公式完整指南，第 211 頁）。但是，當我嘗試使用 MC 方法計算此選項的值時，我得到的值始終不正確。
下面是我的 MC 模擬程式碼。首先是對 SDE 進行數值積分的函式：
# Euler scheme for two GBMs (no correlation) with same drift and volatility # Returns the terminal value (prices at last time step) gbm <- function(mu, sigma, max_time, num_steps, init_value){ h = max_time / num_steps paths <- matrix(NA, num_steps+1, 2) paths[1, ] = init_value normals = matrix(rnorm(num_steps*2, sd=sqrt(h)), num_steps, 2) for (i in 1:num_steps){ paths[i+1, ] = paths[i, ] + (mu * paths[i, ] * h) + (sigma * paths[i, ] * normals[i, ]) } return(paths[num_steps, ]) }然後是蒙地卡羅方法。請注意，我計算了最高看漲期權的價格和普通看漲期權的價格：
trials <- 10000 maxes <- array(NA, trials) max_payoffs <- array(NA, trials) vanilla_payoffs <- array(NA, trials) for(i in 1:trials){ # Compute terminal values of the SDEs terminal_values <- gbm(mu=0.02, sigma=0.2, max_time=3, num_steps=1000, init_value=c(1, 1)) # Vanilla call payoff just on one of the GBM - for assuring my numerical integration correct vanilla_payoffs[i] <- max(terminal_values[1] - 1, 0) # Call on the maximum of the two assets - strike 1 maxes[i] = max(terminal_values) max_payoffs[i] = max(maxes[i] - 1, 0) } # Mean of the payoffs + 95% confidence interval mean(max_payoffs) * exp(-0.02 * 3) sd(max_payoffs * exp(-0.02 * 3)) * 2 / sqrt(trials) # Mean of the vanilla call payoffs mean(vanilla_payoffs) * exp(-0.02 * 3)對於最多兩個資產的呼叫，我的樣本平均值是 $ 0.2839 \pm 0.0064 $ 這與正確的值相差甚遠 $ 0.2235 $ . 但是我的普通看漲期權幾乎完全正確 $ 0.1656 $ 與真實值相比 $ 0.1646 $ .
為了清楚起見，參數是 $ \sigma=0.2 $ , $ r=0.02 $ , $ S^{1,2}_0 = 1 $ , $ K=1 $ , $ T=3 $ , $ \rho=0 $ .
如果有人能解釋我哪裡出錯了，我將不勝感激。
編輯：根據@Yoda And Friends 的回答，我添加了不使用數值積分的python 程式碼。但它仍然給出不正確的價格：
def terminal_spots(trials, r, sigma, t, spot): normals = np.random.normal(size = (trials, 2)) return spot * np.exp(t * (r - 0.5 * sigma * sigma) + sigma * np.sqrt(t) * normals)和
def mc_call_max_two_assets(trials, r, sigma, t, spot, strike): terminals = terminal_spots(trials, r, sigma, t, spot) max_terminal = terminals.max(1) payoffs = np.maximum(max_terminal - strike, 0) mn = payoffs.mean() * np.exp(-r*t) conf_interval = (payoffs * np.exp(-r*t)).std() * 2 / np.sqrt(trials) return mn, conf_interval
我想您要計算以下收益： $$ \pi_T = \max\left[ \max(S_T^1, S_T^2) - K, 0 \right] $$ 如果動態用以下動態表示（從您的程式碼中，應該是這種情況）： $$ dS_t^i = rS_t^idt + \sigma S_t^idW_t^i $$ 在哪裡 $ W = (W^1, W^2) $ 是二維 SBM，那麼，我想你可以稍微簡化一下你的程式碼！
根據您對問題的定義，或有索賠與路徑無關。這意味著您只對查找感興趣 $ S_T^i $ 這只是一個隨機變數。讓我們看看能不能找到它的分佈！進行以下變換： $$ d(\log(S_t^i)) = (r - \frac{1}{2}\sigma^2)dt + \sigma dW_t^i $$ 這只是以下積分錶示的簡寫符號： $$ \log(S_t^i) - \log(S_0^i) = \int_0^t (r - \frac{1}{2}\sigma^2)ds + \int_0^t \sigma dW_s^i $$ 存在 $ r $ 和 $ \sigma $ 常數，並且由於 $ \int_0^t \sigma dW_s^I \sim N(0, \int_0^t \sigma^2) = N(0, \sigma^2 t) $ 我們有： $$ \log(S_t^i) \sim N\left(\log(S_0^i) + (r - \frac{1}{2}\sigma^2)t, \ \ \sigma^2 t\right) $$ 然後我們可以“模擬”一個標準的正態分佈 $ Z \sim N(0, 1) $ 使用我們的筆記型電腦： $$ \log(S_t^i) \sim \log(S_0^i) + (r - \frac{1}{2}\sigma^2)t + \sigma \sqrt{t} Z $$ 考慮到這一點： $ S_t = \exp(\log(S_t)) $ . 所以我們的蒙地卡羅方法歸結為：
def generateStock(initialValue, sigma, rfr, t, randomGenerator): z = RandomGenerator.nextGaussian(0, 1) return initialValue*Math.exp((rfr - 0.5*sigma*sigma)*t + sigma*Math.sqrt(t)*z) def monteCarloMaxOption(initialValue, sigma, rfr, strike, maturity, numberOfSimulations): sample = [0]*numberOfSimulations #Just create an array of length nos for w in range(numberOfSimulations){ stock1 = generateStock(initialValue, sigma, rfr, maturity) stock2 = generateStock(initialValue, sigma, rfr, maturity) maximumStock = Math.max(stock1, stock2) sample[w] = Math.max(maximumStock - strike, 0.0) } discountFactor = Math.exp(-rfr*maturity) priceToday = mean(sample)*discountFactor sampleError = errorFunction(sample) return priceToday, error我很抱歉使用了虛擬碼（受 Java 影響的 Python 風格…：D），但我發現 R 真的很難！
如果回報不是你想要的，請告訴我！！！
我認為您的分析定價公式有問題：
我為 ( p. 211 ) 上指定的分析定價公式提供了一個 R 腳本，它給出的呼叫價格為 0.2853。
library(pbivnorm) maxassets_analytical <- function(r, T, K, sigma1, sigma2, S1, S2, b1, b2, rho){ y1 <- (log(S1/K) + (b1 + sigma1^2/2)*T)/ (sigma1*sqrt(T)) y2 <- (log(S2/K) + (b2 + sigma2^2/2)*T)/ (sigma2*sqrt(T)) sig <- sqrt(sigma1^2 + sigma2^2 - 2 * sigma1 * sigma2 * rho) d <- (log(S1/S2) + (b1 - b2 + sig^2/2)*T)/ (sig*sqrt(T)) rho1 <- (sigma1 - rho * sigma2)/sig rho2 <- (sigma2 - rho * sigma1)/sig call <- S1 * exp((b1 - r) * T) * pbivnorm(y1, d, rho1) + S2 * exp((b2 - r) * T) * pbivnorm(y2, -d + sig * sqrt(T), rho2) - K * exp(-r * T) * (1 - pbivnorm(-y1+sigma1*sqrt(T), -y2 + sigma2 * sqrt(T), rho)) return(list("call price" = call, "y1" = y1, "y2" = y2)) }在哪裡 $ b_1 $ 和 $ b_2 $ 分別是持有成本， $ S_t^1 $ 和 $ S_t^2 $ . 對於不支付股息的股票，持有成本等於無風險利率。插入您指定的參數以及 $ b_1=b_2=0.02 $ 給出0.2853：
為簡潔起見，我提供了以下公式的圖片。
這裡， $ M(a,b,\rho) $ 是累積二元正態分佈函式（見第 13 章）。
確認：
我通過實現相應看跌期權的分析定價公式來驗證我的程式碼，因為他提供了（第 212 頁）中指定的範例計算。不僅如此，看跌期權還嵌套了看漲期權的分析公式，這提供了一種驗證方法：
我實現了 put 定價功能，
R程式碼可以在下面找到。按照他的例子並因此插入參數， $ S_1=100 $ , $ S_2=105 $ , $ X=98 $ , $ T=0.5 $ , $ r=0.05 $ , $ b_1=-0.01 $ , $ b_2=-0.04 $ , $ \sigma_1 = 0.11 $ , $ \sigma_2=0.16 $ , $ \rho = 0.63 $ ，我得到以下結果：
這與他的範例中的結果一致：
我希望這有幫助！
看跌期權程式碼：
maxassets_analyticalput <- function(r, T, K, sigma1, sigma2, S1, S2, b1, b2, rho){ sig <- sqrt(sigma1^2 + sigma2^2 - 2 * sigma1 * sigma2 * rho) d <- (log(S1/S2) + ((b1 - b2) + (sig^2)/2)*T)/ (sig*sqrt(T)) cmax0 <- S2 * exp( (b2 - r) * T) + S1 * exp( (b1 - r) * T) * pnorm(d) - S2 * exp((b2 - r)*T) * pnorm(d-sig*sqrt(T)) cmax <- maxassets_analytical(r, T, K, sigma1, sigma2, S1, S2, b1, b2, rho) pmax <- K * exp(-r * T) - cmax0 + cmax[[1]] return(list("put price" = pmax, "cmax, K = 0" = cmax0, "cmax" = cmax[[1]], "d" = d, "sigma" = sig, "y1" = cmax[[2]], "y2" = cmax[[3]])) }
