Ornstein-Uhlenbeck Process

Ornstein-Uhlenbeck Process
Ornstein-Uhlenbeck Process

Introducción

El proceso de Ornstein–Uhlenbeck fue introducido en 1930 por los físicos Leonard Ornstein y George Eugene Uhlenbeck en el contexto del estudio del movimiento browniano con fricción. A diferencia del proceso de Wiener, que modela trayectorias sin memoria ni tendencia, el proceso de Ornstein–Uhlenbeck incorpora un mecanismo de reversión a la media, lo que lo hace más realista para fenómenos físicos y económicos. En física, describe la velocidad de una partícula sometida a fuerzas aleatorias y fricción en un fluido. En finanzas, se utiliza ampliamente para modelar tasas de interés (modelo de Vasicek), spreads, volatilidad estocástica y precios que tienden a un nivel de equilibrio. En biología, aparece en modelos de evolución de rasgos bajo selección estabilizadora. En ingeniería y machine learning, se usa para generar ruido correlacionado en algoritmos de control y aprendizaje por refuerzo. Su importancia radica en que es un proceso gaussiano, markoviano y estacionario en régimen límite, lo que permite obtener resultados analíticos para media, varianza y autocovarianza. Esto lo convierte en una herramienta fundamental tanto teórica como aplicada en múltiples disciplinas donde existen dinámicas con tendencia al equilibrio.

Definición

El proceso de Ornstein–Uhlenbeck es un proceso estocástico continuo \(\{X_t\}_{t \geq 0}\) definido como la solución de la ecuación diferencial estocástica:

\[ dX_t = \theta(\mu - X_t)\,dt + \sigma\,dW_t \]

donde \(\theta > 0\) es el parámetro de velocidad de reversión a la media, \(\mu \in \mathbb{R}\) es la media de largo plazo, \(\sigma > 0\) es el parámetro de volatilidad, y \(W_t\) es un movimiento browniano estándar.

Este proceso describe una dinámica en la cual la variable \(X_t\) tiende a regresar hacia el nivel \(\mu\) con una fuerza proporcional a la desviación \((\mu - X_t)\), mientras está sujeta a perturbaciones aleatorias modeladas por el término \(\sigma dW_t\).

Demostración

Reescribimos: \[ dX_t + \theta X_t dt = \theta \mu dt + \sigma dW_t \]

Multiplicamos por el factor integrante \(e^{\theta t}\): \[ d(e^{\theta t} X_t) = \theta \mu e^{\theta t} dt + \sigma e^{\theta t} dW_t \]

Integramos de \(0\) a \(t\): \[ e^{\theta t} X_t = X_0 + \theta \mu \int_0^t e^{\theta s} ds + \sigma \int_0^t e^{\theta s} dW_s \]

La solución de la ecuación diferencial estocástica está dada por:

\[ X_t = X_0 e^{-\theta t} + \mu(1 - e^{-\theta t}) + \sigma \int_0^t e^{-\theta (t-s)} dW_s \]

Propiedades principales

Media: \[ \mathbb{E}[X_t] = \mu + (X_0 - \mu)e^{-\theta t} \]

Varianza: \[ \mathbb{V}ar(X_t) = \frac{\sigma^2}{2\theta}(1 - e^{-2\theta t}) \]

Autocovarianza (para \(s \le t\)):

\[ \mathbb{C}ov(X_s, X_t) = \mathrm{Var}(X_s)\, e^{-\theta (t-s)} \]

Y en forma general

\[ \mathbb{C}ov(X_s, X_t) = \frac{\sigma^2}{2\theta} e^{-\theta |t-s|} \]

El proceso de Ornstein–Uhlenbeck es un proceso gaussiano y markoviano con reversión a la media. El parámetro \(\theta\) controla la rapidez con la que el proceso regresa a \(\mu\), mientras que \(\sigma\) determina la magnitud de las fluctuaciones aleatorias. A largo plazo, el proceso alcanza un régimen estacionario con distribución normal:

\[ X_t \sim \mathcal{N}\left(\mu, \frac{\sigma^2}{2\theta}\right) \]

Media

Sabemos que la solución de la ecuación diferencial estocástica está dada por:

\[ X_t = X_0 e^{-\theta t} + \mu(1 - e^{-\theta t}) + \sigma \int_0^t e^{-\theta (t-s)} dW_s \]

Tomando esperanza: \[ \mathbb{E}[X_t] = \mathbb{E}\left[X_0 e^{-\theta t} + \mu(1 - e^{-\theta t}) + \sigma \int_0^t e^{-\theta (t-s)} dW_s \right] \]

por propiedad de la esperanza \[ \mathbb{E}[X_t] = \mathbb{E}\left[X_0 e^{-\theta t}\right] + \mathbb{E}\left[\mu(1 - e^{-\theta t})\right] + \mathbb{E}\left[\sigma \int_0^t e^{-\theta (t-s)} dW_s \right] \]

por propiedad fundamental de Ito \(\mathbb{E}[\int_0^t f(s) \, dW_s] = 0\), dado que el término \(\sigma \int_0^t e^{-\theta (t-s)} dW_s\) depende de la derivada del moviemiento Browniano, entonces, su valor es cero,

\[ \mathbb{E}[X_t] = \mathbb{E}\left[X_0 e^{-\theta t}\right] + \mathbb{E}\left[\mu(1 - e^{-\theta t})\right] \] y, al ser constantes ambos términos, se tiene la expresión análitica del valor esperado del proceso de Ornstein–Uhlenbeck

\[ \boxed{\mathbb{E}[X_t] = X_0 e^{-\theta t} + \mu(1 - e^{-\theta t})} \]

Varianza

La expresión de la varianza es \[ \mathbb{V}ar(X_t) = \mathbb{E}[X_t^2] - (\mathbb{E}[X_t])^2 \]

Para el proceso de Ornstein–Uhlenbeck \[ \mathbb{V}ar(X_t) = \mathbb{E} \left[\left(X_0 e^{-\theta t} + \mu(1 - e^{-\theta t}) + \sigma \int_0^t e^{-\theta (t-s)} dW_s\right)^2 \right] - \\ \left(\mathbb{E}\left[X_0 e^{-\theta t} + \mu(1 - e^{-\theta t}) + \sigma \int_0^t e^{-\theta (t-s)} dW_s \right] \right)^2 \]

No voy a escribir todo el desarrollo de forma explícita, solo comentar que del segundo término ya sabemos el resultado \(\mathbb{E}[X_t] = X_0 e^{-\theta t} + \mu(1 - e^{-\theta t})\), si observamos detenidamente, al realizar los cálculos se van a eliminar los primeros dos términos de la expresión \(\mathbb{E}[X_t^2]\), con los de \(\mathbb{E}[X_t]\), y los términos que dependen de \(\sigma \int_0^t e^{-\theta (t-s)} dW_s\) al obtener el valor esperado van a ser cero, excepto su cuadrado, el único término que sobrevive, es decir

\[ \mathbb{V}ar(X_t) = \mathbb{E} \left[\sigma^2 \left( \int_0^t e^{-\theta (t-s)} dW_s \right)^2 \right] \]

\[ \mathbb{V}ar(X_t) = \sigma^2 \mathbb{E} \left[ \int_0^t \int_0^t e^{-\theta (t-s)} e^{-\theta (t-u)} dW_s dW_u \right] \]

\[ \mathbb{V}ar(X_t) = \sigma^2 \int_0^t \int_0^t e^{-\theta (t-s)} e^{-\theta (t-u)} \mathbb{E} \left[dW_s dW_u \right] \]

Si \(s=u\), entonces \(\mathbb{E} \left[dW_s dW_u \right]= dW_s\), si fuesen distintas el resultado sería cero, entonces

\[ \mathbb{V}ar(X_t) = \sigma^2 \int_0^t e^{-2\theta (t-s)} dW_s \]

Ahora, toca resolver la integral para ello hacemos cambio de variable \(u = t - s\), entonces, \(\quad dW_u = - dW_s\), y

\[ \int_0^t e^{-2\theta (t-s)} \, dW_s = \int_t^0 e^{-2\theta u} \, (-dW_u) = \int_0^t e^{-2\theta u} \, dW_u = \frac{1 - e^{-2\theta t}}{2\theta} \]

por tanto, la varianza del proceso Ornstein–Uhlenbeck es

\[ \boxed{ \mathbb{V}ar(X_t) = \frac{\sigma^2}{2\theta} \left(1 - e^{-2\theta t} \right) } \]

Autocovarianza

\[ \mathbb{C}ov(X_s,X_t) = \mathbb{E}\big[(X_s - \mathbb{E}[X_s])(X_t - \mathbb{E}[X_t])\big] \]

Para \(s \le t\), y recordando el resultado del valor esperado

\[ \mathbb{E}[X_t] = X_0 e^{-\theta t} + \mu(1 - e^{-\theta t}) \]

entonces \[ X_t - \mathbb{E}[X_t] = \sigma \int_0^t e^{-\theta (t-u)} dW_u \] \[ X_s - \mathbb{E}[X_s] = \sigma \int_0^s e^{-\theta (s-v)} dW_v \]

por lo tanto

\[ \mathbb{C}ov(X_s,X_t) = \mathbb{E}\left[\left(\sigma \int_0^s e^{-\theta (t-u)} dW_u \right) \left(\sigma \int_0^t e^{-\theta (s-v)} dW_v \right)\right] \]

\[ \mathbb{C}ov(X_s,X_t) = \mathbb{E}\left[\sigma^2 \int_0^s \int_0^t e^{-\theta (t-u)} e^{-\theta (s-v)} dW_u dW_v \right] \] y \[ \mathbb{C}ov(X_s,X_t) = \sigma^2 \int_0^s \int_0^t e^{-\theta (t-u)} e^{-\theta (s-v)} \mathbb{E}\left[ dW_u dW_v \right] \]

Usando las propiedades de Itô, \(\mathbb{E}\left[ dW_u dW_v \right]=dW_u\) si \(u=v\), y cero si \(u \neq v\), entonces

\[ \mathbb{C}ov(X_s,X_t) = \sigma^2 \int_0^s e^{-\theta (t-u)} e^{-\theta (s-u)} dW_u \]

simplificando el integrando \(e^{-\theta (s-u)} e^{-\theta (t-u)} = e^{-\theta (s+t-2u)}= e^{-\theta (s+t)} e^{-\theta 2u}\), entonces

\[ \mathbb{C}ov(X_s,X_t) = \sigma^2 \int_0^s e^{-\theta (s+t)} e^{-\theta 2u} dW_u \]

\[ \mathbb{C}ov(X_s,X_t) = \sigma^2 e^{-\theta (s+t)} \int_0^s e^{-\theta 2u} dW_u \]

resolviendo la integral \(\int_0^s e^{2\theta u} \, du = \left[ \frac{e^{2\theta u}}{2\theta} \right]_0^s = \frac{e^{2\theta s} - 1}{2\theta}\), tenemos

\[ \mathbb{C}ov(X_s,X_t) = \sigma^2 e^{-\theta (s+t)} \left( \frac{e^{2\theta s} - 1}{2\theta} \right) \]

simpleficando

\[ \boxed{ \mathbb{C}ov(X_s, X_t) = \frac{\sigma^2}{2\theta} \left( e^{-\theta (t-s)} - e^{-\theta (t+s)} \right), \quad s < t } \]

podemos probar que para \(t < s\),

\[ \mathbb{C}ov(X_s, X_t) = \frac{\sigma^2}{2\theta} \left( e^{-\theta (s-t)} - e^{-\theta (t+s)} \right), \quad t < s \]

por tanto, en forma general para \(s>0\), y \(t>0\), se tiene

\[ \boxed{ \mathbb{C}ov(X_s, X_t) = \frac{\sigma^2}{2\theta} \left( e^{-\theta |t-s|} - e^{-\theta (t+s)} \right), \quad s>0, \quad t>0 } \]

Simulación del Proceso de Ornstein-Uhlenbeck

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
from dbconnection import MySQLDatabase
import statsmodels.api as sm
sns.set_theme()

La siguiente función nos va a permitir relizar múltiples trayectorias y con un número definido de epasos, considerando también los parámetros del proceso Ornstein-Uhlenbeck revisado anteriormente, y retornarnos un dataframe.


class StochasticSimulation():

    def __init__(self, T=None, steps=100, n_times=1):
        self.steps = steps
        self.n_times = n_times

        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)

    # proceso Ornstein-Uhlenbeck
    def OrnsteinUhlenbeck(self, X0, theta, mu, sigma):
        ou = np.zeros((self.n_times, self.steps + 1))

        for i in range(self.n_times):
            ou[i, 0] = X0
            for j in range(1, self.steps + 1):
                Z = np.random.normal(0, 1)
                ou[i, j] = (
                    ou[i, j-1] * np.exp(-theta * self.dt)
                    + mu * (1 - np.exp(-theta * self.dt))
                    + sigma * np.sqrt((1 - np.exp(-2 * theta * self.dt)) / (2 * theta)) * Z
                )

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

La función simula una trayectoria del proceso de Ornstein-Uhlenbeck utilizando los parámetros especificados (\(X_0\), \(\theta\), \(\mu\), \(\sigma\)). Calcula la media teórica y la desviación estándar del proceso a lo largo del tiempo para poder compararlas con la simulación.Genera los intervalos de confianza del \(95\%\) basados en la desviación estándar teórica. Finalmente grafica la trayectoria simulada junto con los límites superior e inferior del intervalo de confianza, lo cual permite visualizar la reversión a la media y la dispersión esperada del proceso en función del tiempo.

# parámetros
N = 10000        # número de pasos
M = 1            # número de trayectorias simuladas
X0=-10
theta= 0.0006
mu=-0.02
sigma=0.05

# simulación
sim = StochasticSimulation(steps=N, n_times=M)
df = sim.OrnsteinUhlenbeck(X0,theta,mu,sigma)


# Colormap
cmap = cm.get_cmap("Grays", M)
# lista de t
t = sim.t
T = sim.T
std_teorica = np.sqrt((sigma**2 / (2 * theta)) * (1 - np.exp(-2 * theta * t)))
# Media teórica
mean_B =  X0 * np.exp(-theta * t) + mu * (1 - np.exp(-theta * 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"Proceso  Ornstein-Uhlenbeck [{M} trayectoria]")
plt.xlabel("Tiempo t")
plt.ylabel("S(t)")
plt.grid(True)
plt.legend()
plt.show()

Ahora vemos otra simulación modificando los parámetros de Ornstein-Uhlenbeck.

N = 10000        # número de pasos
M = 1            # número de trayectorias simuladas
X0=10
theta= 0.0009
mu= 2
sigma=0.05

# simulación
sim = StochasticSimulation(steps=N, n_times=M)
df = sim.OrnsteinUhlenbeck(X0,theta,mu,sigma)


# Colormap
cmap = cm.get_cmap("Grays", M)
# lista de t
t = sim.t
T = sim.T
std_teorica = np.sqrt((sigma**2 / (2 * theta)) * (1 - np.exp(-2 * theta * t)))
# Media teórica
mean_B =  X0 * np.exp(-theta * t) + mu * (1 - np.exp(-theta * 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"Proceso  Ornstein-Uhlenbeck [{M} trayectoria]")
plt.xlabel("Tiempo t")
plt.ylabel("S(t)")
plt.grid(True)
plt.legend()
plt.show()

Ahora vamos a simular múltiples trayectorias del proceso de Ornstein–Uhlenbeck utilizando los parámetros especificados (\(X_0\), \(\theta\), \(\mu\), \(\sigma\)) y un número determinado de pasos (\(N\)) y trayectorias (\(M\)). El gráfico muestra las trayectorias simuladas junto con los límites superior e inferior del intervalo de confianza, mostrando la dispersión y variabilidad de las trayectorias, lo cual permite visualizar la reversión a la media, la dispersión esperada y el comportamiento colectivo del proceso a través de múltiples simulaciones.

# parámetros
N = 10000        # número de pasos
M = 1000           # número de trayectorias simuladas
X0=-10
theta= 0.0006
mu=-0.02
sigma=0.05

# simulación
sim = StochasticSimulation(steps=N, n_times=M)
df = sim.OrnsteinUhlenbeck(X0,theta,mu,sigma)

# gráfico
cmap = cm.get_cmap("Grays", M*2)
# lista de t
t = sim.t
T = sim.T
std_teorica = np.sqrt((sigma**2 / (2 * theta)) * (1 - np.exp(-2 * theta * t)))
# Media teórica
mean_B =  X0 * np.exp(-theta * t) + mu * (1 - np.exp(-theta * 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"Proceso  Ornstein-Uhlenbeck [{M} trayectoria]")
plt.xlabel("Tiempo t")
plt.ylabel("S(t)")
plt.grid(True)
plt.legend()
plt.show()

Vamos seleccionar ciertos pasos específicos (\(10\), \(100\), \(1000\), y \(10000\)) de la simulación del proceso y construiimos histogramas de las posiciones correspondientes. Los gráficos permiten visualizar cómo evolucionan las distribuciones a lo largo del tiempo, y apreciar la reversión a la media. En pasos pequeños (por ejemplo \(10\) o \(100\)) los histogramas muestran que las posiciones están muy concentradas alrededor del valor inicial, todavía no se ha alcanzado el equilibrio. En pasos grandes (\(1000\) o \(10000\)) los histogramas tienden a una distribución más amplia y estable, que refleja la distribución estacionaria teórica del proceso.

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

Veamos el comportamiento para la segunda gráfica mostrada, pero ahora con múltiples trajectorias. Al igual que el caso anterior podemos observar la reversión a la media, la dispersión esperada y el comportamiento colectivo del proceso a través de múltiples simulaciones.

# parámetros
N = 10000        # número de pasos
M = 1000          # número de trayectorias simuladas
X0=10
theta= 0.0009
mu= 2
sigma=0.05

# simulación
sim = StochasticSimulation(steps=N, n_times=M)
df = sim.OrnsteinUhlenbeck(X0,theta,mu,sigma)

# gráfico
cmap = cm.get_cmap("Grays", M*2)
# lista de t
t = sim.t
T = sim.T
std_teorica = np.sqrt((sigma**2 / (2 * theta)) * (1 - np.exp(-2 * theta * t)))
# Media teórica
mean_B =  X0 * np.exp(-theta * t) + mu * (1 - np.exp(-theta * 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"Proceso  Ornstein-Uhlenbeck[{M} trayectoria]")
plt.xlabel("Tiempo t")
plt.ylabel("S(t)")
plt.grid(True)
plt.legend()
plt.show()

Vamos a visualizar cómo evolucionan las distribuciones a lo largo del tiempo, y apreciar la reversión a la media. A medida que el tiempo avanza, la distribución se centra alrededor de la media de largo plazo \((\mu)\), evidenciando que el proceso tiende a regresar a su valor medio. Por otro lado, la dispersión de las posiciones aumenta con el tiempo hasta estabilizarse, mostrando cómo la varianza crece hacia un valor límite dado por la teoría.

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

Proceso Ornstein-Uhlenbeck en dos dimensiones

Ahora vamos a graficar el Ornstein-Uhlenbeck, no sobre el parámetro \(t\), sino visualizarlo en 2 dosmensiones.Veamos el caso de una trayectoria en 2 dimensiones.

# parámetros
N = 10000        # número de pasos
M = 1000           # número de trayectorias simuladas
X0=-10
theta= 0.0006
mu=-0.02
sigma=0.05

# simulación
sim = StochasticSimulation(steps=N, n_times=M)
df = sim.OrnsteinUhlenbeck(X0,theta,mu,sigma)

# Colormap (gradiente de gris)
cmap = cm.get_cmap("Greys", M*2)
centro = (0, 0)
circulo = patches.Circle(centro, 20, 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 < 1:
        ax.plot(df.iloc[i][:], df.iloc[i + 1][:], linewidth=0.2, color=cmap(i))
    else:
        continue
plt.title('Simulación de un proceso  Ornstein-Uhlenbeck: {} steps and {} paths'.format(N, M))
plt.xlabel('X(t)')
plt.ylabel('X(t)')
plt.gca().set_xticklabels([])  # quitar etiquetas eje X
plt.gca().set_yticklabels([])  # quitar etiquetas eje Y
plt.grid(True)
plt.show()

Ahora el mismo caso, pero ahora agregando \(800\) trayectorias

# parámetros
N = 10000        # número de pasos
M = 1000           # número de trayectorias simuladas
X0=0
theta= 0.0006
mu= 0
sigma=0.09

# simulación
sim = StochasticSimulation(steps=N, n_times=M)
df = sim.OrnsteinUhlenbeck(X0,theta,mu,sigma)

# Colormap (gradiente de gris)
cmap = cm.get_cmap("Greys", M*2)
centro = (0, 0)
circulo = patches.Circle(centro, 15, 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 de un proceso  Ornstein-Uhlenbeck: {} steps and {} paths'.format(N, M))
plt.xlabel('X(t)')
plt.ylabel('x(t)')
plt.gca().set_xticklabels([])
plt.gca().set_yticklabels([])
plt.grid(True)
plt.show()

Estimación de parámetros

La estimación de los parámetros del proceso de Ornstein-Uhlenbeck se realiza a partir de una serie temporal observada \(\{X_t\}\), como log-precios o spreads financieros. El proceso Ornstein-Uhlenbeck se define mediante la ecuación diferencial estocástica: \[ dX_t = \theta (\mu - X_t) \, dt + \sigma \, dW_t, \] cuya solución

\[ X_t = X_0 e^{-\theta t} + \mu(1 - e^{-\theta t}) + \sigma \int_0^t e^{-\theta (t-s)} dW_s \]

donde \(\theta\) es la velocidad de reversión a la media, \(\mu\) la media de largo plazo y \(\sigma\) la volatilidad del proceso.

Para estimar los parámetros, se puede discretizar la ecuación como un modelo AR(1): \[ X_{t+\Delta t} = \alpha + \beta X_t + \varepsilon_t, \] con \(\varepsilon_t \sim \mathcal{N}(0, \sigma_t^2)\). Los coeficientes \(\alpha=\mu(1 - e^{-\theta t})= \mu(1 - \beta)\) y \(\beta=e^{-\theta t}\) se obtienen mediante regresión lineal ordinaria.

A partir de esta regresión, los parámetros continuos se calculan como: \[ \theta = -\frac{\ln \beta}{\Delta t}, \quad \mu = \frac{\alpha}{1 - \beta}, \quad \sigma = \sigma_t \sqrt{\frac{2\theta}{1 - \beta^2}}. \]

Esta metodología permite calibrar el proceso de Ornstein-Uhlenbeck, simular trayectorias y construir intervalos de confianza teóricos \[ \mathbb{E}[X_t] = X_0 e^{-\theta t} + \mu (1 - e^{-\theta t}), \quad \mathbb{V}ar(X_t) = \frac{\sigma^2}{2 \theta} \left( 1 - e^{-2 \theta t} \right), \] lo que facilita el análisis de la reversión a la media y la dispersión esperada del proceso.

Aplicación a la serie MSTF

Veamos un ejemplo práctivco, vamoa a descargar series de precios díarios de la compañia MSTF de enero de \(2010\) a marzo de \(2026\).

# conección a base de datos
db = MySQLDatabase("financialmarkets")
query = """  
    SELECT sp.date,
           sp.company_id, 
           c.symbol,
           sp.close_price
    FROM stock_prices_sp500 as sp
    LEFT JOIN companies as c
        ON sp.company_id = c.company_id
    WHERE DATE(date) >= '2010-01-01' AND symbol = 'MSFT';
"""
# ejecuta consulta
base = db.execute_query(query)
## ✅ Conexión exitosa
base.shape
## (5922, 4)

Ahora veamos el gráfico de la serie. Se observa una tendencia general al alza, con crecimiento gradual hasta \(2019\) y un aumento más pronunciado entre \(2020\) y \(2025\). También se aprecian picos y caídas abruptas, lo que indica volatilidad temporal significativa. Hacia el final de la serie, el precio presenta una corrección o retroceso, evidenciando fluctuaciones propias del mercado. En general, la serie refleja tendencia, volatilidad y posibles episodios de reversión a la media, lo cual justifica la aplicación del proceso de Ornstein-Uhlenbeck.

df = base[['date', 'close_price']]
df['date'] = pd.to_datetime(df['date'])
df = df.set_index('date')
log_price = np.log(df.astype(float))
# gráfico
plt.figure(figsize=(14,6))
plt.plot(log_price.index, log_price['close_price'])
plt.title('Precio de cierre logarítmico de MSFT')
plt.xlabel('Fecha')
plt.ylabel('Log Precio')
plt.show()

De acuerdo a lo comentado en al inicio de la sección, vamos estimar los parámetros del proceso Ornstein-Uhlenbeck. La siguientes líneas de código realiza la estimación de los parámetros de un proceso de a partir de una serie temporal de log_price. Primero, se crean los pares \((X_t, X_{t+1})\) necesarios para construir un modelo AR(1) discreto y se eliminan los valores faltantes en caso de existir. Luego ajusta una regresión lineal ordinaria de \(X_{t+1}\) en función de \(X_t\) con constante, obteniendo los coeficientes \(\alpha\) y \(\beta\). A partir de estos coeficientes, calcula los parámetros continuos del proceso Ornstein-Uhlenbeck: la velocidad de reversión a la media (\(\theta\)), la media de largo plazo (\(\mu\)) y la volatilidad (\(\sigma\)), usando la varianza de los residuos del modelo.

# Cargar serie
X = log_price  

# Crear pares (Xt, Xt+1)
df['X_t'] = X
df['X_t1'] = X.shift(-1)

# elimina valores na
df_model = df.dropna()

# modelo
X_t_const = sm.add_constant(df_model['X_t'])
model = sm.OLS(df_model['X_t1'], X_t_const).fit()
# parámetros
alpha = model.params[0]
beta = model.params[1]

# Parámetros OU
dt = 1  # diario: 1/252
theta = -np.log(beta) / dt
mu = alpha / (1 - beta)

# calculo de Sigma
residuals = model.resid
var_eps = np.var(residuals)
sigma = np.sqrt(var_eps * (2 * theta) / (1 - np.exp(-2 * theta * dt)))

# Resultados
print(f"theta: {round(theta,2)} \n mu: {round(mu,2)} \n sigma: {round(sigma,2)}")
## theta: 0.0 
##  mu: 16.29 
##  sigma: 0.01

Finalemnet, vamos a graficar el ajuste con la serie real. El gráfico muestra un ajuste de un proceso de Ornstein-Uhlenbeck a la serie de precios de la acción MSFT.La línea negra representa la serie real del logaritmo de precios de MSFT, y la línea azul muestra la serie ajustada teóricamente por el proceso Ornstein-Uhlenbeck. La banda azul clara alrededor de la serie ajustada representa el intervalo de confianza teórico del \(95\%\), esto indica la zona en la que el modelo predice que la serie debería fluctuar con alta probabilidad. La serie real tiende a oscilar dentro del intervalo de confianza, lo que sugiere que el proceso Ornstein-Uhlenbeck captura razonablemente bien la dinámica general de la acción

# parámetros
N = df.shape[0]        # número de pasos
M = 1                  # número de trayectorias simuladas
# redondeamos valor inicial de la serie 
X0= round(log_price['close_price'].to_list()[0], 2)
# valores estimados del paso anterior
theta= theta
mu= mu
sigma=sigma

# simulación
sim = StochasticSimulation(steps=N, n_times=M)
df = sim.OrnsteinUhlenbeck(X0,theta,mu,sigma)

# Colormap
cmap = cm.get_cmap("Grays", M)
# lista de t
t = sim.t
T = sim.T
std_teorica = np.sqrt((sigma**2 / (2 * theta)) * (1 - np.exp(-2 * theta * t)))
# Media teórica
mean_B =  X0 * np.exp(-theta * t) + mu * (1 - np.exp(-theta * 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))
plt.plot(log_price['close_price'].reset_index(drop=True), label='Serie real')
for i in range(M):
    plt.plot(t, df.iloc[i,:].to_list(), color=cmap(i), linewidth=0.5, label='Trayectoria simulada')
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"Ajuste de un Proceso  Ornstein-Uhlenbeck a la serie MSFT")
plt.xlabel("Tiempo t")
plt.ylabel("S(t)")
plt.grid(True)
plt.legend()
plt.show()

Conclusiones

El Proceso de Ornstein-Uhlenbeck es un ejemplo clásico de proceso estocástico estacionario de media-reversión, ampliamente utilizado para modelar fenómenos donde las variables tienden a regresar hacia un valor promedio a lo largo del tiempo. Finalmente, la aplicación práctica demuestra que el proceso Ornstein-Uhlenbeck, aunque simple, captura aspectos esenciales de series financieras o de tiempo que exhiben tendencia a la media, ofreciendo un marco robusto para análisis cuantitativo y simulación de escenarios. La combinación de teoría, simulación y ajuste empírico garantiza que las decisiones basadas en este modelo tengan un fundamento sólido tanto matemático como práctico.

Referencias

  • Wolfgang, Paul and Jörg Baschnagel(2013). Stochastic Processes: From Physics to Finance. 2da edition. Springer.
  • 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.
  • Lefebvre Mario(2007). Applied Stochastics Processes. Springer, New York.