Geometric Brownian Motion

Movimiento Browniano Geométrico

Definición:

Un proceso estoc'astico \(S(t)\) se dice que sigue un Movimiento Browniano Geométrico(MBG) si satisface la siguiente ecuación diferencial estocástica:

\[\begin{equation} dS(t) = \mu S(t)dt + \sigma S(t) dW(t), \ \ S(0) = S_0 \end{equation} \]

donde \(W(t)\) es un proceso de Wiener, \(\mu\) y \(\sigma\) son constantes.

Poposición (lema de Itô)

Sea X(t) un proceso estocástico que satisface la siguiente ecuaci'on diferencial estocástica

\[\begin{eqnarray*} dX(t) = a(X(t), t) dt + b(X(t), t)dW(t) \end{eqnarray*} \]

y sea \(f(x,t)\) una funci'on doblemente diferenciaciable con variables \(x\) y \(t\) el lema de It^o establece que

\[ \begin{eqnarray*} df(X(t), t) = \left[\frac{\partial f(X,t)}{\partial t} + a(X,t)\frac{\partial f(X,t)}{\partial x} + \frac{b^2(X,t)}{2} \frac{\partial^2 f(X,t)}{\partial x^2}\right]dt + b(X,t) \frac{\partial f(X,t)}{\partial x} dW(t) \end{eqnarray*} \]

Aplicando el Lema de Itô y solucionando, se tiene \[ \begin{eqnarray*} \frac{dS(t)}{S(t)} & = & \mu dt + \sigma dW(t)\\ d lnS(t) & = & \mu dt + \sigma dW(t)\\ d lnS(t) & = & \mu dt + \sigma dW(t) - \frac{1}{2} \sigma^2 dt\\ d lnS(t) & = & \left(\mu - \frac{1}{2} \sigma^2 \right) dt + \sigma dW(t)\\ \int_0^T d lnS(t) & = & \int_0^T \left(\mu - \frac{1}{2} \sigma^2 \right) dt + \int_0^T \sigma dW(t)\\ ln S(T) - ln S(0) & = & \left(\mu - \frac{\sigma^2}{2} \right)T + \sigma W(T)\\ ln S(T) & = & ln S(0) + \left(\mu - \frac{\sigma^2}{2} \right)T + \sigma W(T) \end{eqnarray*} \]

donde \[ \begin{eqnarray*} S(T) & = & S_0 \exp \left(\left(\mu - \frac{\sigma^2}{2} \right)T + \sigma W(T) \right)\\ \end{eqnarray*} \]

Proposición (Distribución lognormal)

Una variable aleatoria \(X\) sigue una distribución lognormal, con \(media = \mu\) y \(varianza = \sigma^2\), si \(X = e^Y\), donde \(Y\) es una variable aleatoria con distribuci'on normal. Y tiene la siguinete función de distribución: \[ \begin{eqnarray} f(x) & = \frac{1}{x\sqrt{2\pi \sigma}} \exp \left(\frac{1}{2\sigma^2} \left(ln(x)-\mu \right)^2 \right), & x>0 \\ & = 0, & d.o.m \nonumber \end{eqnarray} \]

Demostración:

Supongamos que \(Y \sim N(\mu,\sigma)\), entonces empleando el método de la función de distribución se tiene que: \[ \begin{eqnarray*} p(e^y \leq x) & = & p(y \leq ln(x))\\ & = & F_y(ln(x))\\ & = & f_y(ln(x))\frac{1}{x}\\ & = & \frac{1}{x\sqrt{2\pi \sigma}} \exp \left(\frac{1}{2\sigma^2} \left(ln(x)-\mu \right)^2 \right), \qquad x> 0 \end{eqnarray*} \]

Derivado del resultado anterior y de las propiedades de la distribución normal, la distribución de \(S(t)\) viene dada por

\[ \begin{align*} \ln S(t) &= \ln S(0) + \left( \mu - \frac{\sigma^2}{2} \right)t + \sigma W(t) \sim N\left( \ln S(0) + \left( \mu - \frac{\sigma^2}{2} \right)t, \, \sigma^2 t \right), \\ S(t) &= \exp\!\left\{ \ln S(0) + \left( \mu - \frac{\sigma^2}{2} \right)t + \sigma W(t) \right\} \sim LN\left( \ln S(0) + \left( \mu - \frac{\sigma^2}{2} \right)t, \, \sigma^2 t \right) \end{align*} \]

Ya que conocemos la distribución de \(S(t)\) emplearemos el siguiente resultado para hallar la media y varianza del proceso.

Media y varianza del MBG

Si un proceso estoc'astico \(S(t)\) sigue un Movimiento Browniano Geom'etrico(MBG) (descrito anteriormente), entonces, su esperanza y varianza vienen dadas por: \[ \begin{eqnarray*} \mathbb{E}[X] & = & S_0\exp \left(\mu t \right), y\\ \mathbb{V}ar(X) & = & S_0^2 \exp \left( 2\mu t \right) \left( \exp \left( \sigma^2 t \right) -1 \right) \end{eqnarray*} \]

Demostración:

Esperanza de la distribución Lognormal

La esperanza viene definida por \[ \begin{align*} \mathbb{E}[X] &= \int_{0}^{\infty} \frac{x}{x \sqrt{2\pi \sigma^2}} e^{-\frac{1}{2\sigma^2}(\ln x - \mu)^2} \, dx \\[8pt] \text{Sea } w = \frac{\ln x - \mu}{\sigma}, \text{ donde } x = e^{\sigma w + \mu} \text{ y } dx = x\sigma \, dw. \\[8pt] \mathbb{E}[X] &= \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi \sigma^2}} e^{-\frac{1}{2}w^2} \, x\sigma \, dw \\[6pt] &= \frac{\sigma}{\sqrt{2\pi \sigma^2}} \int_{-\infty}^{\infty} e^{-\frac{1}{2} w^2} e^{\sigma w + \mu} \, dw \\[6pt] &= \int_{-\infty}^{\infty} \frac{e^{\mu + \frac{\sigma^2}{2}}}{\sqrt{2\pi}} e^{-\frac{1}{2}(w - \sigma)^2} \, dw \\[6pt] &= e^{\mu + \frac{\sigma^2}{2}} \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi \sigma^2}} e^{-\frac{1}{2}(w - \sigma)^2} \, dw \end{align*} \]

dado que \(\int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2} \left(w-\sigma \right)^2} dw \sim N(0,1)\) e integra \(1\), entonces \[ \begin{eqnarray*} \mathbb{E}[X] = e^{\mu + \frac{\sigma^2}{2}} \end{eqnarray*} \]

Varianza de la distribución Lognormal

El segundo momento viene definido como: \[ \begin{align*} \mathbb{E}[X^2] &= \int_{0}^{\infty} \frac{x^2}{x \sqrt{2\pi \sigma^2}} e^{-\frac{1}{2\sigma^2}(\ln x - \mu)^2} \, dx \\[8pt] \text{realizando el mismo cambio de variable, se tiene que:} \\[6pt] \mathbb{E}[X^2] &= \int_{-\infty}^{\infty} \frac{e^{2\mu + 2\sigma^2}}{\sqrt{2\pi \sigma^2}} e^{-\frac{1}{2}w^2 + 4w\sigma - 4\sigma^2} \, dw \\[6pt] &= \int_{-\infty}^{\infty} \frac{e^{2\mu + 2\sigma^2}}{\sqrt{2\pi \sigma^2}} e^{-\frac{1}{2}(w - 2\sigma)^2} \, dw \\[6pt] &= e^{2\mu + 2\sigma^2} \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi \sigma^2}} e^{-\frac{1}{2}(w - 2\sigma)^2} \, dw \end{align*} \]

dado que \(\int_{-\infty}^{\infty} e^{-\frac{1}{2} \left(w -2\sigma \right)^2} dw \sim N(2\sigma,1)\) e integra \(1\), y

\[ \begin{eqnarray*} \mathbb{E}[x^2] & = & e^{2\mu + 2\sigma^2} \end{eqnarray*} \]

entonces, la varianza viene dada por \[ \begin{align*} \mathbb{V}ar(X) &= E[X^2] - \left(E[X]\right)^2 \\[6pt] &= e^{2\mu + 2\sigma^2} - \left(e^{\mu + \frac{1}{2}\sigma^2}\right)^2 \end{align*} \]

Estimación por máxima verosimilitud

Sea \(Xi = logS(t_i) - logS(t_{i-1})\). Para encontrar los parámetros \(\mu\) y \(\sigma\), primero hallamos función de verosimilitud \[ \begin{eqnarray*} L(x_i) & = & \prod_{i=1}^{n} \frac{1}{x_i\sqrt{2\pi \sigma^2 t}} \exp \left(\frac{1}{2\sigma^2t} \left(x_i- \left(\mu - \frac{1}{2} \sigma^2 \right)t\right)^2 \right)\\ & = & \left[\frac{1}{x_i 2\pi \sigma^2 t} \right]^{n/2} \exp \left(\frac{1}{2\sigma^2t} \sum_{i=1}^{n} \left(x_i - \left(\mu - \frac{1}{2} \sigma^2 \right) t\right)^2 \right) \end{eqnarray*} \] la log-verosimilitud viene dada por \[ \begin{eqnarray*} l(\mu,\sigma) & = & - \frac{n}{2} ln(x_i 2\pi \sigma^2 t) - \frac{1}{2\sigma^2 t} \sum_{i=1}^{n} \left(x_i- \left(\mu - \frac{1}{2} \sigma^2 \right) t\right)^2 \end{eqnarray*} \]

sea \(v= \left(\mu - \frac{1}{2} \sigma^2 \right) t\) y \(w=\sigma^2t\), entonces obteniendo las derivadas con respecto a los par'ametros \[ \begin{eqnarray*} \frac{\partial l}{\partial v} & = & \frac{1}{w} \sum_{i=1}^{n} \left(x_i- v \right)\\ \frac{\partial l}{\partial w} & = & - \frac{n}{2w} + \frac{1}{2w^2} \sum_{i=1}^{n} \left(x_i- v \right)^2 \end{eqnarray*} \] Igualando a cero se obtiene \[ \begin{eqnarray*} \sum_{i=1}^{n} x_i- n\hat{v} & = & 0\\ - \frac{n}{2\hat{w}} + \frac{1}{2\hat{w}^2} \sum_{i=1}^{n} \left(x_i- \hat{v} \right)^2 & = & 0 \end{eqnarray*} \]

solucionando \[ \begin{eqnarray*} \hat{v} & = & \frac{1}{n}\sum_{i=1}^{n}x_i, y\\ \hat{w} & = & \frac{1}{n}\sum_{i=1}^{n} \left(x_i - \hat{v}\right)^2 \end{eqnarray*} \] dado que definimos anteriormente \(v= \left(\mu - \frac{1}{2} \sigma^2 \right) t\) y \(w=\sigma^2t\), entonces \[ \begin{eqnarray*} \hat{\mu} & = & \frac{1}{\hat{\sigma}^2} + \frac{1}{nt} \hat{v}, y\\ \hat{\sigma^2} & = & \frac{1}{nt} \hat{w} \end{eqnarray*} \]

Además de la estimación por máxima verosimilitud, exiten otras técnicas para poder estimar el valor de los parámetros, por ejemplo, modelos Bayesianos, Por otro lado, se han aplicado modelos de deeplearning para resolver ecuaciones diferenciales estocásticas.

Simulación

Este código define una clase \(StochasticSimulation()\) para simular procesos estocásticos en tiempo discreto. El método init inicializa el número de pasos, repeticiones (\(n_times\)), el tiempo total \(T\) y genera un vector de tiempos \(t\) con incremento \(dt\). El método \(GeometricBrownianMotion(S0,alpha,beta)\) el cual simula trayectorias. Se realizan \(n_times\) simulaciones en paralelo, almacenando cada trayectoria en filas del arreglo. Finalmente, convierte el resultado en un DataFrame de pandas.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import matplotlib.cm as cm
import matplotlib.patches as patches
sns.set_theme()

class StochasticSimulation():

    def __init__(self, T=None, steps=100, n_times=1):
        self.steps = steps
        self.n_times = n_times
        # Si T no se proporciona, usar T = steps
        if T is None:
            self.T = steps
        else:
            self.T = T

        self.dt = self.T/ steps
        self.t = np.linspace(0, self.T, steps+1)

    # Movimiento Browniano Geométrico
    def GeometricBrownianMotion(self, S0, alpha, beta):
        gbm = np.zeros((self.n_times, self.steps + 1))

        for i in range(self.n_times):      
            gbm[i, 0] = S0
            for j in range(1, self.steps + 1):
                B = np.random.normal(0, self.dt)
                gbm[i, j] = gbm[i, j-1] * np.exp((alpha - 0.5 * beta**2) * self.dt + beta * B)

        df = pd.DataFrame(gbm, columns=np.round(self.t, 3))
        return df

En este ejemplo, primero se definen los parámetros de la simulación: se tomarán \(𝑁=10,000\) pasos, lo que permitirá generar una trayectoria muy detallada del proceso, y se simulará \(𝑀=1\), una única trayectoria.

A continuación, se crea una instancia de la clase \(StochasticSimulation()\), pasando como argumentos el número de pasos \(N\)y el número de trayectorias \(M\) a simular. Esta clase se encarga de manejar los cálculos necesarios para generar procesos estocásticos.

Finalmente, se llama al método \(GeometricBrownianMotion(S0,mu,sigma)\) de la instancia, recibe como argumento \(S0\) valor inicial, asó como \(alpha\) y \(beta\), con los cuales nos genera la trayectoria simulada. El resultado se almacena en df, un DataFrame que contiene los valores de la trayectoria a lo largo de los pasos de tiempo definidos.

Simulación de un Movimiento Browniano Geométrico en una dimensión

Simulación de una trajectoria

# Parámetros
# Parámetros
N = 10000        # número de pasos
M = 1           # número de trayectorias simuladas
S0=10
mu=0.0001
sigma=0.003

sim = StochasticSimulation(steps=N, n_times=M)
df = sim.GeometricBrownianMotion(S0, alpha=mu, beta=sigma)

# Colormap
cmap = cm.get_cmap("Grays", M*2)
# lista de t
t = sim.t
T = sim.T
std_teorica = S0*np.sqrt(np.exp(2*mu*t)*(np.exp(sigma**2*t) -1))
# Media teórica
mean_B = S0*np.exp(mu*t)
# Intervalos de confianza 95%
conf_upper =  mean_B + 1.96 * std_teorica
conf_lower =  mean_B - 1.96 * std_teorica
# Gráfico
plt.figure(figsize=(14,6))
for i in range(M):
    plt.plot(t, df.iloc[i,:].to_list(), color=cmap(i), linewidth=0.5)
plt.plot(t, conf_upper, 'blue', linewidth=0.3, linestyle='--')
plt.plot(t, conf_lower, 'blue', linewidth=0.3, linestyle='--')
plt.fill_between(t, conf_lower, conf_upper, color='lightblue', alpha=0.2, label='Intervalo de confianza 95% (teórico)')
plt.title(f"Movimiento Browniano Geométrico[{M} trayectoria]")
plt.xlabel("Tiempo t")
plt.ylabel("S(t)")
plt.grid(True)
plt.legend()
plt.show()

Simulación de múltiples trajectorias

# Parámetros
N = 10000        # número de pasos
M = 1000           # número de trayectorias simuladas
S0=10
mu=0.0001
sigma=0.003

sim = StochasticSimulation(steps=N, n_times=M)
df = sim.GeometricBrownianMotion(S0, alpha=mu, beta=sigma)

# Colormap
cmap = cm.get_cmap("Grays", M*2)
# lista de t
t = sim.t
T = sim.T
std_teorica = S0*np.sqrt(np.exp(2*mu*t)*(np.exp(sigma**2*t) -1))
# Media teórica
mean_B = S0*np.exp(mu*t)
# Intervalos de confianza 95%
conf_upper =  mean_B + 1.96 * std_teorica
conf_lower =  mean_B - 1.96 * std_teorica
# Gráfico
plt.figure(figsize=(14,6))
for i in range(M):
    plt.plot(t, df.iloc[i,:].to_list(), color=cmap(i), linewidth=0.5)
plt.plot(t, conf_upper, 'blue', linewidth=0.3, linestyle='--')
plt.plot(t, conf_lower, 'blue', linewidth=0.3, linestyle='--')
plt.fill_between(t, conf_lower, conf_upper, color='lightblue', alpha=0.2, label='Intervalo de confianza 95% (teórico)')
plt.title(f"Movimiento Browniano Geométrico[{M} trayectoria]")
plt.xlabel("Tiempo t")
plt.ylabel("S(t)")
plt.grid(True)
plt.legend()
plt.show()

Veamos las distribuciones a distintos tiempos

# Tiempos a graficar
steps_to_plot = [10, 100, 1000, 10000]

# Crear figura con 2x2 subplots
fig, axes = plt.subplots(2, 2, figsize=(10, 6))
axes = axes.flatten()  # para indexar fácilmente

for i, step in enumerate(steps_to_plot):
    data = df[step]
    mean = data.mean()
    var = data.var()

    sns.histplot(data, bins=30, kde=True, ax=axes[i])
    axes[i].set_title(f'Distribución {step}')
    axes[i].set_xlabel('Posición')
    axes[i].set_ylabel('Frecuencia')
    axes[i].legend([f'Media = {mean:.2f}, Varianza = {var:.2f}'])

plt.tight_layout()
plt.show()

Simulación de un Movimiento Browniano Geométrico en dos dimensiones

Finalmente, tomamos los datos de la simulación anterior y gráficamos pares de caminatas aleatorias para ver su evolución en dos dimensiones.

# Colormap (gradiente de gris)
cmap = cm.get_cmap("Greys", M*2)
centro = (30, 30)
circulo = patches.Circle(centro, 40, edgecolor='grey', facecolor='black', alpha=0.8, lw=1, ls='--')

# Crear la figura y los ejes
fig, ax = plt.subplots(figsize=(14, 10))
ax.add_patch(circulo)
# Brownian Motion
for i in range(df.shape[0] - 1):
    if i < 800:
        ax.plot(df.iloc[i][:], df.iloc[i + 1][:], linewidth=0.2, color=cmap(i))
    else:
        continue
plt.title('Simulación Movimiento Broaniano Geométrico: {} steps and {} paths'.format(N, M))
plt.xlabel('S(t)')
plt.ylabel('S(t)')
plt.gca().set_xticklabels([])  # quitar etiquetas eje X
plt.gca().set_yticklabels([])  # quitar etiquetas eje Y
plt.grid(True)
plt.show()

Referencias

Ubbo F Wiersema (2008). Brownian Motion Calculus. John Wiley & Sons Ltd.
Sheldon M. Ross(1995). Stochastic Processes. Wiley Series in Probability and Statistics: Probability and Statistics. John Wiley & Sons.
LefebvreMario(2007). Applied Stochastics Processes. Springer, New York.