如何為期貨價差期權定價？

假設我有兩個期貨合約 $ F_1(0,T) $ 和 $ F_2(0,T) $ 在兩個不同的相關基礎上。

如果我假設兩個標的都遵循具有波動性的 GBM $ \sigma_1 $ 和 $ \sigma_2 $ 分別，並且他們的回報是相關的 $ \rho $ , 是否有一個封閉式公式可以用來對價差的歐式期權進行定價 $ X_{1,2} = \max( F_1(S,T) - F_2(S,T), 0 ) $ 在某個到期日 $ S \leq T $ ?

否則，我們通常如何對這些東西進行定價？蒙地卡羅使用風險中性動態和貼現收益？

我所知道的最好的近似值是 Li, Deng and Zhou, 2006。這是一種解析近似值，其中價格以直接公式表示，因此易於實現。

如果你想非常準確，這是我的論文 J Choi (2018) ( Arxiv )。它處理任何資產線性組合的期權，例如一籃子期權和亞洲期權以及點差期權。有一些關於其他解析近似的性能的討論。

正如 Jaehyuk Choi 所說，他有一個非常準確的點差期權 + 籃子 + 亞洲期權定價器。他沒有提到的是他已經在這里為你用 R 編碼了它：https ://github.com/PyFE/SumBSM-R/blob/master/SumBSM/blksmd.R

他指出 Li, Deng, Zhou 2006 是一個準確且易於編碼的程式碼。他已經對它進行了編碼（對於看漲期權），我只是添加了看跌期權平價來計算看跌期權價值（請注意，它會給出帶有負數的錯誤結果，也許有人可以編寫實際公式，因為我無法訪問論文）。你通過點喜歡：spot <- c(spot1, spot2)和卷喜歡vol <- c(vol1, vol2)- 我根本不是 R 使用者，但仍然可以稍微遵循程式碼。程式碼沒有處理一些邊界“中斷”條件，但至少這可以讓您領先一步：

blksmd_CalcSpreadOptLDZ <- function( strk, spot, t.exp, vol, rho, callput, r = 0, d = 0 ){
 call <- rep(NA, length(strk))
 fwd <- spot * exp((r-d)*t.exp)
 std <- vol * sqrt(t.exp)
 
 rho.comp <- sqrt(1-rho*rho)*std[1]
 mu1 <- log(fwd[1]) - 0.5*std[1]*std[1]
 mu2.exp <- fwd[2] * exp(-0.5*std[2]*std[2]) # R = exp(mu2) with y0 = 0

 # vector for the rest
 r.plus.k <- mu2.exp + strk # (R+K) in DLZ (vectorized)

 epsilon <- -1/(2*rho.comp) * (std[2]*std[2]*mu2.exp*strk)/(r.plus.k*r.plus.k)
 
 C3 <- ( mu1 - log(r.plus.k) )/rho.comp
 D3 <- ( rho*std[1] - std[2]*mu2.exp/r.plus.k )/rho.comp
 
 C2 <- C3 + std[2]*( D3 + epsilon*std[2] )
 D2 <- D3 + 2*std[2]*epsilon
 
 C1 <- C3 + rho*std[1]*(D3 + epsilon*rho*std[1]) + rho.comp
 D1 <- D3 + 2*rho*std[1]*epsilon

 I_S1 <- blksmd_CalcSpreadOptLDZ_I( C1, D1, epsilon )
 I_S2 <- blksmd_CalcSpreadOptLDZ_I( C2, D2, epsilon )
 I_K  <- blksmd_CalcSpreadOptLDZ_I( C3, D3, epsilon )

 price.fwd <- fwd[1]*I_S1 - fwd[2]*I_S2 - strk*I_K

 if (callput) {
     return(exp(-r*t.exp)*price.fwd)
     } else {
     return((exp(-r*t.exp)*price.fwd) + exp(-r*t.exp)*strk - (spot[1]-spot[2]))
     }  
}

#' An auxilary function used by blksmd_CalcSpreadOptLDZ
#' 
blksmd_CalcSpreadOptLDZ_I <- function( u, v, eps ) {
 u2 <- u*u
 u4 <- u2*u2
 v2 <- v*v
 v4 <- v2*v2
 v6 <- v4*v2
 v2.sqrt <- sqrt(1+v2)
 
 arg.norm <- u/v2.sqrt
 J0 <- pnorm(arg.norm)
 J1 <- (1+(1+u2)*v2)/v2.sqrt^5L * dnorm(arg.norm)
 J2 <- (6*(1-u2)*v2 + (21-2*u2-u4)*v4 + 4*(3+u2)*v6 - 3)*u/v2.sqrt^11L * dnorm(arg.norm)
 
 return( J0 + eps*(J1 + 0.5*eps*J2))
}

所以這是他的程式碼（添加了一個 puts），如果你喜歡，你應該支持他的答案。在上面的連結中，他還編寫了非常複雜的程序，謝謝 Jaehyuk ！！！

引用自：https://quant.stackexchange.com/questions/15961