Movimiento browniano geométrico

Los precios cambian en porcentaje

El movimiento browniano aritmético añade el mismo ruido en dólares cueste la acción 5 o 500, y puede bajar de cero. Una acción se comporta distinto: un movimiento del 1 % es un 1 % a cualquier precio. El movimiento browniano geométrico pone el azar en el cambio porcentual: dS = mu * S * dt + sigma * S * dW. Es el modelo de Black–Scholes. mu es el rendimiento esperado anual y sigma la volatilidad, la magnitud que se cotiza en el 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

La solución exacta y el menos un medio sigma al cuadrado

Resolver la ecuación da S(t) = S(0) * exp((mu - sigma**2 / 2) * t + sigma * W(t)). ¿Por qué ese extraño - sigma**2 / 2? Porque la exponencial es convexa: por la desigualdad de Jensen, el ruido eleva el promedio de exp(ruido), y restar sigma**2 / 2 lo cancela exactamente, de modo que el precio *esperado* crece a tasa mu. Las *trayectorias* típicas crecen más despacio, a mu - sigma**2 / 2, por eso el precio mediano está por debajo de la media. El siguiente tema, el lema de Itô, deriva este término.

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

Muchas trayectorias: la media y la mediana

Simula miles de trayectorias de un año y mira dónde terminan. El precio final *medio* se acerca a S0 * exp(mu * T), la fórmula del precio esperado. La *mediana* es menor, S0 * exp((mu - sigma**2/2) * T), y con alta volatilidad la diferencia es grande: la mayoría de las trayectorias acaba por debajo de la media, impulsada por unas pocas espectaculares. Esta asimetría a la derecha, una distribución lognormal, es la forma de los precios de las acciones en horizontes largos.

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 y sigma

A partir de precios muestreados cada dt años, los rendimientos logarítmicos r = diff(log(prices)) son normales independientes de media (mu - sigma**2/2) * dt y desviación típica sigma * sqrt(dt). Así sigma = std(r) / sqrt(dt) y mu = mean(r) / dt + sigma**2 / 2. La estimación de la volatilidad es buena con unos cientos de puntos; la deriva es notoriamente ruidosa, porque su error estándar solo disminuye con el *tiempo* total, no con el número de observaciones.

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)); los cambios porcentuales son aleatorios y los precios siguen siendo positivos.
• El - sigma**2 / 2 anula la convexidad de exp: el precio esperado es S0 * exp(mu * t), pero la trayectoria mediana crece a mu - sigma**2 / 2.
• Estime con los rendimientos log: sigma = std(r) / sqrt(dt), mu = mean(r) / dt + sigma**2 / 2; la volatilidad se estima bien, la deriva mal.

A continuación: el lema de Itô, la regla de la cadena que explica de dónde sale el - sigma**2 / 2.

Lab: qué se califica

• Dar un paso exacto con gbm_step25 pts
• Simular una trayectoria de precios con gbm_path25 pts
• Promediar todos los futuros con gbm_expected25 pts
• Recuperar los parámetros con estimate_gbm25 pts