Como usar a simulação de Monte Carlo com movimento browniano

Neste artigo implementamos de forma simples uma simulação de Monte Carlo de preços de ações com movimento browniano geométrico (GBM), explicando a lógica por trás para que você possa executá-la corretamente.

Em finanças matemáticas, um ativo costuma ser modelado como um movimento browniano geométrico, ou seja, seu preço obedece a uma equação diferencial estocástica (EDE). O GBM supõe que o preço sobe e desce de forma aleatória, mas com uma tendência geral de alta ou baixa ao longo do tempo. Ele captura isso com dois ingredientes: a velocidade média de variação do preço (a deriva) e o quanto ele oscila aleatoriamente (a volatilidade). Também supõe que os retornos logarítmicos do ativo seguem uma distribuição normal.

Neste exemplo usamos o preço da ação da Apple e executamos uma simulação de Monte Carlo com GBM no mesmo período. Em seguida, plotamos as simulações junto com a trajetória real da ação.

Duas formas de modelar o GBM

1. A abordagem de Black-Scholes

O modelo de Black-Scholes supõe que os retornos de um ativo seguem uma distribuição log-normal. Aplicando a matemática, chega-se à fórmula padrão do preço do ativo:

S(t) = S(0) * exp((μ - 0.5 * σ²) * t + σ * W(t))

Aqui, μ é a deriva (o retorno médio anualizado da ação), σ é a volatilidade (o desvio padrão dos retornos) e W(t) é um movimento browniano padrão. O termo -0.5 * σ² é uma correção para o efeito de “arrasto” nos retornos, e já faz parte da fórmula.

Simulação de Monte Carlo com GBM para a ação da Apple (2021-2024)

Vamos à prática. O código abaixo simula o preço da Apple de 2021 a 2024 em Python.

import numpy as np
import pandas as pd
import yfinance as yf
import matplotlib.pyplot as plt

apple = yf.download("AAPL", start="2021-01-01", end="2024-01-01")
apple.columns = apple.columns.levels[0]
price = apple.Close
returns = price.pct_change().dropna()

t = len(returns)           # number of daily steps
dt = 1                     # one trading day per step
mu = returns.mean()        # drift
sigma = returns.std()      # volatility
simulations = 1000
S0 = price.iloc[0]

def gbm_paths(mu, sigma, W):
    return np.exp((mu - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * W)

Ws = np.random.standard_normal(size=(simulations, t))
gbm = gbm_paths(mu, sigma, Ws)
gbm = np.insert(gbm, 0, S0, axis=1).cumprod(axis=1)

x_axis = range(0, t + 1)
fig, ax = plt.subplots(figsize=(12, 5.5))
for path in gbm:
    ax.plot(x_axis, path, color="grey", linewidth=0.1)
ax.plot(x_axis, price, color="red", linewidth=1)

2. A abordagem de transformação logarítmica

Se os retornos são log-normais, basta tomar o logaritmo dos retornos. Depois dessa transformação, os dados podem ser modelados com um movimento browniano padrão. Chega-se ao mesmo modelo de fundo, visto por uma lente matemática diferente.

log_returns = np.log(returns + 1)
mu_log = log_returns.mean()
sigma_log = log_returns.std()

def bm_paths(mu, sigma, W):
    return mu + sigma * W

bm = bm_paths(mu_log, sigma_log, Ws)
bm = S0 * np.exp(bm.cumsum(axis=1))
bm = np.insert(bm, 0, S0, axis=1)

fig, ax = plt.subplots(figsize=(12, 5.5))
for path in bm:
    ax.plot(x_axis, path, color="grey", linewidth=0.1)
ax.plot(x_axis, price, color="red", linewidth=1)

Comparando os dois resultados

Os dois gráficos são praticamente idênticos. Com retornos logarítmicos perfeitamente normais, eles convergem para o mesmo resultado, o que confirma que ambos os métodos são formas válidas de executar uma simulação de Monte Carlo com movimento browniano geométrico.