Responder com médias de sorteios aleatórios
Muitas perguntas financeiras não têm fórmula arrumada: a chance de uma carteira perder 20% em um ano, o preço de uma opção cujo pagamento depende da trajetória inteira. Monte Carlo as responde na força bruta. Escreva a quantidade como uma média, E[f(X)], simule muitos sorteios de X, aplique f a cada um e tire a média. Pela lei dos grandes números, a média amostral converge para a esperança verdadeira. Uma probabilidade é apenas a esperança de um indicador sim/não, então "a fração de simulações em que o evento ocorreu" a 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 preços e probabilidades
Para simular uma ação, sorteie choques normais padrão e passe-os pela fórmula lognormal do tópico anterior: os n preços terminais saem de uma única expressão sobre arrays. Depois qualquer pergunta sobre a distribuição dos preços vira uma média sobre o array. A probabilidade de terminar abaixo de 80 é np.mean(prices < 80), o preço esperado é prices.mean() e um quantil das perdas é np.quantile(...). Passe o gerador rng para a função em vez de criá-lo dentro, para que quem chama controle a repetibilidade.
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)Precificar uma opção
Uma call europeia paga max(S - K, 0) no vencimento. No mundo neutro ao risco todo ativo cresce à taxa livre de risco rate, então o preço justo é o pagamento esperado descontado: simule preços terminais com mu = rate, calcule o pagamento de cada um, tire a média e multiplique por exp(-rate * t). Para essa opção simples existe uma fórmula fechada (Black-Scholes, 10,4506 para os números abaixo) e a simulação chega perto. O ponto do Monte Carlo é que as mesmas poucas linhas ainda funcionam para pagamentos sem nenhuma 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: escreva a quantidade como uma média
E[f(X)], simule muitos sorteios, apliquef, tire a média. - Uma probabilidade é a média de uma comparação:
np.mean(prices < 80); passe o gerador como argumento para que a execução seja repetível. - Preço da opção =
exp(-rate * t) * mean(max(S - K, 0))com os preços simulados amu = rate; informe também o erro.
A seguir: o erro de Monte Carlo e a convergência, quão precisa é uma estimativa e como obter mais com menos.
Lab: o que é avaliado
- Atirar dardos num quadrado para achar pi com
estimate_pi25 pts - Contar a frequência de um evento com
estimate_probability25 pts - Simular onde uma ação poderia terminar com
simulate_terminal_prices25 pts - Precificar uma opção por simulação com
mc_call_price25 pts