使用 R 的蒙地卡羅期權定價
我正在嘗試使用 R 使用 Monte Carlo 實現一個普通的歐洲期權定價器。下面是我的程式碼，用於在股票遵循 GBM 的情況下對非股息支付股票的歐洲普通普通看漲期權進行定價。
出於教學原因，我使用了解析公式和 Euler-Maruyama 近似。
但是，將得到的結果與 B&S 模型的結果進行比較，我發現差異很大，因此我想請問您是否可以發現我的 Monte Carlo 程式碼中的錯誤：
# Compute the Black-Scholes European option price on non-dividend paying stock # Setting the B&S parameters value S <- 52 #stock price at time t K <- 50 #strike price tau <- 0.25 #time to maturity T - t (in years) #0.25 = 3 months r <- 0.05 #risk-free annual interest rate sigma <- 0.3 #annual volatility of the stock price (standard deviation) #call B&S fair value d1 <- (log(S/K) + (r + 0.5*sigma^2)*tau)/(sigma*sqrt(tau)) d2 <- d1 - sigma*sqrt(tau) V_BS_Call <- S*pnorm(d1) - K*exp(-r*(tau))*pnorm(d2) #fair value call # Compute the Monte Carlo European option price on non-dividend paying stock # Assuming the non- dividend paying stock follows a Geometric Brownian Motion (GBM) set.seed(2503) #set the seed # Setting the Monte Carlo simulation and GBM parameters tau <- tau #time to expiry (we have already defined this variable) N <- 250 #number of sub intervals dt <- tau/N #length of each time sub interval time <- seq(from=0, to=tau, by=dt) #time moments in which we simulate the process length(time) #it should be N+1 nSim <- 10000 #number of simulations (paths) r <- r #GBM parameter 1 sigma <- sigma #GBM parameter 2 X0 <- S #initial condition (price of the underlying today) #Monte Carlo with analytic formula Z <- matrix(rnorm(nSim*N, mean=0, sd=1),nrow = nSim, ncol = N) #standard normal sample of N elements dW <- Z*sqrt(dt) #Brownian motion increments (N increments)x nSim simulations W <- matrix(numeric(nSim*(N+1)), nrow = nSim, ncol = (N+1)) X_analytic <- numeric(nSim) for(k in 1:nSim){ W[k,] <- c(0, cumsum(dW[k,])) X_analytic[k] <- X0*exp((r - 0.5*sigma^2)*tau + sigma*W[k,ncol(W)]) #Analytic solution } payoff_expiry_call <-pmax(X_analytic-K,0) #pmax preserve the dimension of the matrix, so apply the max function to each element expected_payoff_call <- sum(payoff_expiry_call)/length(payoff_expiry_call) Monte_Carlo_call_price <- exp(-r*(tau))*expected_payoff_call #Monte Carlo with Euler-Maruyama scheme X_EM <- matrix(numeric(nSim*(N+1)), nrow = nSim, ncol = (N+1)) X_EM[,1] <- X0 #first element of X_EM is X0. with the for loop we find the other N elements for(k in 1:nSim){ for(i in 2:ncol(X_EM)){ X_EM[k,i] <- X_EM[k,i-1] + r*X_EM[k,i-1]*dt + sigma*X_EM[k,i-1]*dW[k,i-1] } } payoff_expiry_call <-pmax(X_EM[,ncol(X_EM)]-K,0) #pmax preserve the dimension of the matrix, so apply the max function to each element expected_payoff_call <- sum(payoff_expiry_call)/length(payoff_expiry_call) Monte_Carlo_call_price <- exp(-r*(tau))*expected_payoff_call因此，使用 10,000 次模擬：
- 具有解析公式的蒙地卡羅價格約為 4.535
- 使用 Euler-Maruyama 的 Monte Carlo 價格約為 4.536
- B&S 價格為 4.519
我認為差異太大，但我無法發現錯誤。
您的程式碼看起來不錯，令人鼓舞的是，兩個 MC 模擬都產生了相似的結果。請查看此簡化程式碼，用於蒙地卡羅模擬的分析部分。如你所知，$$ S_T=S_0\exp\left(\left(r-\frac{1}{2}\sigma^2\right)T+\sigma W_T\right). $$呼叫與路徑無關，因此無需模擬整個路徑。我猜你想教你的學生盡可能高效地編碼。自從 $ W_T\sim N(0,T) $ ，可以直接模擬最終的布朗運動。
Z <- rnorm(nSim, mean=0, sd=1) WT <- sqrt(tau) * Z ST = X0*exp((r - 0.5*sigma^2)*tau + sigma*WT) simulated_call_payoffs <- exp(-r*tau)*pmax(ST-K,0) Call_price_MC_anal <- mean(simulated_call_payoffs)如果您稍微玩一下，您確實會獲得與 Black Scholes 封閉形式解決方案不太接近的各種價格。10,000 個樣本值太少，無法準確估計期權價格。嘗試一百萬次模擬。
一般來說，您可以將此作為減少變異數對蒙地卡羅模擬如此重要的動機。該估計值可能是一致且無偏的，但如果您有較大的標準誤差，這對您沒有幫助。回想一下，MC 估計量的信賴區間由下式給出 $$ \hat{C}n \pm z{\delta/2}\frac{s_C}{\sqrt{n}}, $$ 在哪裡 $ \hat{C}_n $ 是估計的贖回價格 $ n $ 模擬和 $ s_c $ 是模擬呼叫值的樣本變異數。顯然，越大 $ n $ ，這個區間越小。如果
nSim=1000000，我得到一個區間 $ [4.51,4.53] $ （BS價格為 $ 4.52 $ ) 但nSim=10000只給出 $ [4.45, 4.69] $ . 95% 信賴區間通過以下方式計算lower_bound <- Call_price_MC_anal - 1.96*sd(simulated_call_payoffs)/sqrt(nSim) upper_bound <- Call_price_MC_anal + 1.96*sd(simulated_call_payoffs)/sqrt(nSim)