Movimento browniano geométrico

Os preços mudam em percentual

O movimento browniano aritmético soma o mesmo ruído em dólares custe a ação 5 ou 500, e pode ficar abaixo de zero. Uma ação se comporta de outro modo: um movimento de 1% é 1% em qualquer preço. O movimento browniano geométrico põe o acaso na variação percentual: dS = mu * S * dt + sigma * S * dW. É o modelo de Black–Scholes. mu é o retorno esperado anual e sigma a volatilidade, a grandeza cotada no mercado.

import math

def gbm_step(s, mu, sigma, dt, z):
    return s * math.exp((mu - 0.5 * sigma ** 2) * dt + sigma * math.sqrt(dt) * z)

print(round(gbm_step(100.0, 0.08, 0.2, 1 / 252, 0.0), 4))     # a calm day: only the (mu - sigma^2/2) drift acts
print(round(gbm_step(100.0, 0.08, 0.2, 1 / 252, 1.0), 4))     # a one-sigma up day
print(round(gbm_step(100.0, 0.08, 0.2, 1 / 252, -50.0), 4))   # an absurd fall: still positive

A solução exata e o menos meio sigma ao quadrado

Resolver a equação dá S(t) = S(0) * exp((mu - sigma**2 / 2) * t + sigma * W(t)). Por que esse estranho - sigma**2 / 2? Porque a exponencial é convexa: pela desigualdade de Jensen, o ruído eleva a média de exp(ruído), e subtrair sigma**2 / 2 o cancela exatamente, de modo que o preço *esperado* cresce à taxa mu. As *trajetórias* típicas crescem mais devagar, a mu - sigma**2 / 2, por isso o preço mediano fica abaixo da média. O próximo tópico, o lema de Itô, deriva esse termo.

import numpy as np

def gbm_path(s0, mu, sigma, dt, z):
    log_moves = (mu - 0.5 * sigma ** 2) * dt + sigma * np.sqrt(dt) * z
    return s0 * np.exp(np.concatenate([[0.0], np.cumsum(log_moves)]))

rng = np.random.default_rng(1)
z = rng.standard_normal(252)
path = gbm_path(100.0, 0.08, 0.2, 1 / 252, z)
print(path[:4].round(3), round(float(path[-1]), 3))

Muitas trajetórias: a média e a mediana

Simule milhares de trajetórias de um ano e veja onde terminam. O preço final *médio* fica perto de S0 * exp(mu * T), a fórmula do preço esperado. A *mediana* é menor, S0 * exp((mu - sigma**2/2) * T), e com volatilidade alta a diferença é grande: a maioria das trajetórias termina abaixo da média, puxada por poucas espetaculares. Essa assimetria à direita, uma distribuição lognormal, é a forma dos preços de ações em horizontes longos.

import numpy as np

rng = np.random.default_rng(5)
s0, mu, sigma, T, steps, paths = 100.0, 0.08, 0.4, 1.0, 252, 20000
dt = T / steps
z = rng.standard_normal((paths, steps))
final = s0 * np.exp(np.sum((mu - 0.5 * sigma ** 2) * dt + sigma * np.sqrt(dt) * z, axis=1))

print(round(float(final.mean()), 1), round(s0 * np.exp(mu * T), 1))                      # close to each other: the expected price
print(round(float(np.median(final)), 1), round(s0 * np.exp((mu - 0.5 * sigma ** 2) * T), 1))   # the median is lower
print(round(float(np.mean(final < s0)), 3))                                              # the share of paths that lose money

Estimar mu e sigma

De preços amostrados a cada dt anos, os retornos logarítmicos r = diff(log(prices)) são normais independentes de média (mu - sigma**2/2) * dt e desvio-padrão sigma * sqrt(dt). Assim sigma = std(r) / sqrt(dt) e mu = mean(r) / dt + sigma**2 / 2. A estimativa da volatilidade é boa com algumas centenas de pontos; a deriva é notoriamente ruidosa, porque seu erro padrão só diminui com o *tempo* total, não com o número de observações.

import numpy as np

def estimate_gbm(prices, dt):
    r = np.diff(np.log(prices))
    sigma = r.std(ddof=1) / np.sqrt(dt)
    return r.mean() / dt + 0.5 * sigma ** 2, sigma

rng = np.random.default_rng(8)
dt = 1 / 252
for years in (1, 10, 100):
    n = int(years * 252)
    log_moves = (0.10 - 0.5 * 0.25 ** 2) * dt + 0.25 * np.sqrt(dt) * rng.standard_normal(n)
    prices = 100 * np.exp(np.concatenate([[0.0], np.cumsum(log_moves)]))
    mu, sigma = estimate_gbm(prices, dt)
    print(years, "years:", round(float(mu), 3), round(float(sigma), 3))           # true mu 0.10 and sigma 0.25

• MBG: S(t) = S(0) * exp((mu - sigma**2 / 2) * t + sigma * W(t)); as variações percentuais são aleatórias e os preços ficam positivos.
• O - sigma**2 / 2 cancela a convexidade de exp: o preço esperado é S0 * exp(mu * t), mas a trajetória mediana cresce a mu - sigma**2 / 2.
• Estime pelos retornos log: sigma = std(r) / sqrt(dt), mu = mean(r) / dt + sigma**2 / 2; a volatilidade se estima bem, a deriva mal.

A seguir: o lema de Itô, a regra da cadeia que explica de onde vem o - sigma**2 / 2.

Lab: o que é avaliado

• Dar um passo exato com gbm_step25 pts
• Simular uma trajetória de preços com gbm_path25 pts
• Fazer a média de todos os futuros com gbm_expected25 pts
• Recuperar os parâmetros com estimate_gbm25 pts