Monte Carlo opción de fijación de precios con R

Question

Estoy tratando de implementar una vainilla opción Europea encargado del precio de Monte Carlo mediante R. En la siguiente no es mi código para la fijación de precios de la unión Europea plain vanilla call opción en el incumplimiento de las acciones que pagan dividendos, bajo el supuesto de que la población sigue una GBM.

Para la enseñanza de razones, he utilizado tanto la fórmula analítica y la de Euler-Maruyama aproximación.

Sin embargo, al comparar los resultados obtenidos con los de B&S modelo, me encontré con una muy gran diferencia, por lo tanto, me gustaría preguntarle si usted puede detectar el error en mi código de 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

Así, utilizando una de 10.000 simulaciones:

• el Monte Carlo precio con fórmula analítica es de aproximadamente 4.535

• el Monte Carlo precio el uso de Euler-Maruyama es de aproximadamente 4.536

• el B&S precio es 4.519

Creo que la diferencia es demasiado grande, pero no puedo descubrir el error.

Answer 1

El código se ve bien y es alentador el hecho de que tanto las simulaciones MC dan resultados similares. Por favor, mire esta simplificado código para la parte analítica de la simulación de Monte Carlo. Como ustedes saben, $$S_T=S_0\exp\left(\left(r-\frac{1}{2}\sigma^2\derecho)T+\sigma W_T\derecho).$$ Una llamada ruta de acceso independiente, así que no hay necesidad para simular la ruta de acceso completa. Supongo que quieres enseñar a tus alumnos a código de la manera más eficiente posible. Desde $W_T\sim N(0,T)$, usted directamente puede simular el final del movimiento Browniano.

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)

Si juegas con un poco de todo, usted realmente obtener diversos precios que no están demasiado cerca de la de Black Scholes forma cerrada de la solución. De 10.000 valores de la muestra son demasiado pocos para estimar con precisión el precio de la opción. Trate de un millón de simulaciones lugar.

Usted podría, en general, el uso de esta una motivación por qué de reducción de varianza es tan crucial para simulaciones de Monte Carlo. La estimación puede ser consistente e imparcial, pero que no ayuda, si usted tiene grandes errores estándar. Recordemos que el intervalo de confianza para la MC estimador está dada por $$ \hat{C}_n \pm z_{\delta/2}\frac{s_C}{\sqrt{n}},$$ donde $\hat{C}_n$ es el estimado de costo de una llamada con $n$ simulaciones y $s_c$ es la varianza de la muestra de la simulación de la llamada de los valores. Obviamente, la mayor $n$, el más pequeño de este intervalo. Si nSim=1000000, puedo obtener un intervalo de $[4.51,4.53]$ (BS precio es de $4.52$) pero nSim=10000 sólo da $[4.45, 4.69]$. El 95% de intervalo de confianza se calcula a través de

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)