Responder con promedios de extracciones aleatorias
Muchas preguntas financieras no tienen una fórmula ordenada: la probabilidad de que una cartera pierda un 20 % en un año, el precio de una opción cuyo pago depende de toda la trayectoria. Monte Carlo las responde por fuerza bruta. Escribe la cantidad como un promedio, E[f(X)], simula muchas extracciones de X, aplica f a cada una y toma la media. Por la ley de los grandes números la media muestral converge a la esperanza verdadera. Una probabilidad es solo la esperanza de un indicador sí/no, así que «la fracción de simulaciones en que ocurrió el evento» la estima.
import math
import numpy as np
# Estimate pi: the quarter circle x^2 + y^2 <= 1 covers pi/4 of the unit square
rng = np.random.default_rng(42)
x, y = rng.random(100000), rng.random(100000)
inside = x ** 2 + y ** 2 <= 1 # True/False for every point
print(round(4 * inside.mean(), 4)) # 3.1378 close to 3.1416, with 100,000 pointsSimular precios y probabilidades
Para simular una acción, extrae choques normales estándar y pásalos por la fórmula lognormal del tema anterior: los n precios terminales salen de una sola expresión sobre arrays. Luego cualquier pregunta sobre la distribución de precios se convierte en un promedio sobre el array. La probabilidad de acabar por debajo de 80 es np.mean(prices < 80), el precio esperado es prices.mean() y un cuantil de las pérdidas es np.quantile(...). Pasa el generador rng a la función en lugar de crearlo dentro, para que quien llama controle la repetibilidad.
def simulate_terminal_prices(rng, s0, mu, sigma, t, n):
z = rng.standard_normal(n)
return s0 * np.exp((mu - 0.5 * sigma ** 2) * t + sigma * math.sqrt(t) * z)
prices = simulate_terminal_prices(np.random.default_rng(1), 100, 0.05, 0.2, 1, 200000)
print(round(float(np.mean(prices < 80)), 4)) # 0.1027 P(price below 80 after a year), about 10%
print(round(float(np.mean(prices > 130)), 4)) # 0.1216 P(price above 130)Valorar una opción
Un call europeo paga max(S - K, 0) al vencimiento. En el mundo neutral al riesgo todo activo crece a la tasa libre de riesgo rate, así que el precio justo es el pago esperado descontado: simula precios terminales con mu = rate, calcula el pago de cada uno, promedia y multiplica por exp(-rate * t). Para esta opción simple existe una fórmula cerrada (Black-Scholes, 10,4506 con los números de abajo) y la simulación se acerca a ella. La gracia de Monte Carlo es que las mismas pocas líneas siguen valiendo para pagos sin ninguna fórmula.
def mc_call_price(rng, s0, strike, rate, sigma, t, n):
prices = simulate_terminal_prices(rng, s0, rate, sigma, t, n) # risk-neutral: the price grows at the risk-free rate
payoff = np.maximum(prices - strike, 0) # max(S - K, 0) for every simulated path
return math.exp(-rate * t) * payoff.mean() # discounted average payoff
price = mc_call_price(np.random.default_rng(1), 100, 100, 0.05, 0.2, 1, 200000)
print(round(price, 4)) # 10.3993 the exact Black-Scholes price is 10.4506: the difference is Monte Carlo error# More points, less error (1 / sqrt(n)): the estimate of pi with the same seed
for n in (100, 10000, 1000000):
g = np.random.default_rng(42)
x, y = g.random(n), g.random(n)
print(n, round(4 * float((x * x + y * y <= 1).mean()), 4))
# 100 3.36
# 10000 3.1488
# 1000000 3.1424- Monte Carlo: escribe la cantidad como un promedio
E[f(X)], simula muchas extracciones, aplicaf, toma la media. - Una probabilidad es la media de una comparación:
np.mean(prices < 80); pasa el generador como argumento para que la ejecución sea repetible. - Precio de la opción =
exp(-rate * t) * mean(max(S - K, 0))con los precios simulados amu = rate; informa también el error.
A continuación: el error de Monte Carlo y la convergencia, qué tan precisa es una estimación y cómo obtener más con menos.
Lab: qué se califica
- Lanzar dardos a un cuadrado para hallar pi con
estimate_pi25 pts - Contar con qué frecuencia ocurre un evento con
estimate_probability25 pts - Simular dónde podría acabar una acción con
simulate_terminal_prices25 pts - Valorar una opción por simulación con
mc_call_price25 pts