程式

兩種資產最大值的看漲期權的蒙地卡羅近似

  • September 9, 2021

我想用收益計算期權的價格 $$ \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 $ ,我得到以下結果:

結果2

這與他的範例中的結果一致

豪格蘇爾

我希望這有幫助!


看跌期權程式碼:

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]]))
}

引用自:https://quant.stackexchange.com/questions/64256