Aproximación de Monte Carlo de la opción de compra sobre el máximo de dos activos

Question

Quiero calcular el precio de la opción con pago \begin{equation} \max \big\{\max\{S^1_T, S^2_T\} - K, 0\big\}, \end{equation} donde $S^{1,2}$ tienen la misma dinámica con 0 correlación. Por lo tanto, \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} donde $W^1$ y $W^2$ son procesos de Wiener independientes bajo la medida de precios $Q$ . Esta opción cuenta con una fórmula analítica de fijación de precios (por ejemplo, en Guía completa de fórmulas de valoración de opciones , p.211). Sin embargo, cuando intento calcular el valor de esta opción utilizando el método MC, obtengo valores que son sistemáticamente incorrectos.

A continuación, mi código para la simulación de MC. Primero una función para hacer la integración numérica de las SDEs:

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

A continuación, el método de Montecarlo. Obsérvese que calculo el precio de una opción de compra sobre el máximo Y el precio de sólo una opción de compra vainilla:

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)

Para la llamada al máximo de dos activos, mi media muestral es $0.2839 \pm 0.0064$ que está muy lejos del valor correcto de $0.2235$ . Sin embargo, mi opción de compra vainilla es casi exactamente correcta $0.1656$ en comparación con el valor real $0.1646$ .

Para que quede claro, los parámetros son $\sigma=0.2$ , $r=0.02$ , $S^{1,2}_0 = 1$ , $K=1$ , $T=3$ , $\rho=0$ .

Estaría muy agradecido si alguien pudiera explicar en qué me estoy equivocando.

EDITAR: He añadido código python que no utiliza integración numérica según la respuesta de @Yoda And Friends. Sin embargo, sigue dando un precio incorrecto:

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)

y

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

Answer 1

Imagino que quieres calcular el siguiente resultado: $$\pi_T = \max\left[ \max(S_T^1, S_T^2) - K, 0 \right]$$ Si la dinámica se expresa con la siguiente dinámica (a partir de su código, debería ser el caso): $$dS_t^i = rS_t^idt + \sigma S_t^idW_t^i$$ donde $W = (W^1, W^2)$ es un SBM bidimensional, entonces, ¡supongo que podrías simplificar un poco tu código!

A partir de su definición del problema, la demanda contingente es independiente del camino. Esto significa que sólo te interesa encontrar $S_T^i$ que es una simple variable aleatoria. Veamos si podemos encontrar su distribución. Tomemos la siguiente transformación: $$d(\log(S_t^i)) = (r - \frac{1}{2}\sigma^2)dt + \sigma dW_t^i$$ que no es más que una notación abreviada para la siguiente representación integral: $$\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$$ En $r$ y $\sigma$ constante, y como $\int_0^t \sigma dW_s^I \sim N(0, \int_0^t \sigma^2) = N(0, \sigma^2 t)$ lo tenemos: $$\log(S_t^i) \sim N\left(\log(S_0^i) + (r - \frac{1}{2}\sigma^2)t, \ \ \sigma^2 t\right)$$ Podemos entonces "simular" una distribución normal estándar $Z \sim N(0, 1)$ utilizando nuestro portátil: $$\log(S_t^i) \sim \log(S_0^i) + (r - \frac{1}{2}\sigma^2)t + \sigma \sqrt{t} Z$$ Teniendo esto en cuenta: $S_t = \exp(\log(S_t))$ . De modo que nuestro método de Monte Carlo se reduce a esto:

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

Siento haber utilizado pseudocódigo (estilo Python con influencias de Java... :D) ¡pero R me parece realmente difícil!

¡¡¡Si el resultado no es lo que pretendía, por favor hágamelo saber!!!

Answer 2

Creo que algo falla en su fórmula de análisis de precios:

He proporcionado un R script para la fórmula analítica de precios especificada en ( p. 211 ) y da un precio de compra de 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))

}

Dónde $b_1$ y $b_2$ es el coste de transporte, respectivamente, $S_t^1$ y $S_t^2$ . En el caso de los valores que no pagan dividendos, el coste del carry es igual al tipo sin riesgo. Introduciendo los parámetros especificados, así como $b_1=b_2=0.02$ da 0,2853:

En aras de la brevedad, he proporcionado una imagen de la fórmula a continuación .

Aquí, $M(a,b,\rho)$ es la función de distribución normal bivariada acumulativa (como se ve en Capítulo 13 ).

---

Verificación:

He verificado mi código aplicando la fórmula analítica de fijación de precios para la opción de venta correspondiente, ya que proporciona un ejemplo de cálculo especificado en ( p. 212 ). No sólo eso, sino que la opción de venta anida la fórmula analítica de la opción de compra, lo que da una forma de verificarla :

He implementado la función de poner precio en R y el código se puede encontrar a continuación. Siguiendo su ejemplo y así insertar los parámetros, $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$ Obtengo los siguientes resultados:

Lo que concuerda con el resultado de su ejemplo :

Espero que esto ayude.

---

Código para la opción de venta:

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