Ecuaciones de Langevin y fluctuaciones alrededor de un equilibrio

Año 2 · Módulo 14

131.1. Relajación y fluctuaciones alrededor de un equilibrio

Cerca de un equilibrio estable, las perturbaciones aleatorias desplazan el estado mientras una fuerza restauradora lo aproxima de nuevo a su valor de referencia. El proceso de Ornstein–Uhlenbeck es una descripción mínima de esa competencia. Su interés está en que tiene una solución probabilística exacta, con la que se calculan las distribuciones transitorias, la varianza estacionaria y la correlación temporal.

Sin fluctuaciones, la ecuación lineal

\frac{dx}{dt}=-\kappa x,\qquad \kappa>0,

tiene la solución $x(t)=x(0)e^{-\kappa t}$ . El parámetro $\kappa$ es una tasa de recuperación y su inverso fija una escala característica de tiempo. El equilibrio $x=0$ es asintóticamente estable, de modo que una perturbación inicial suficientemente representada por este modelo disminuye sin cambiar espontáneamente de dirección. Esa previsibilidad depende de que la ecuación incluya todas las influencias relevantes sobre la variable.

Consideremos una esfera microscópica suspendida en agua. La fricción tiende a reducir su velocidad, mientras las colisiones con numerosas moléculas modifican continuamente su movimiento. Desde la perspectiva de la esfera, las fuerzas microscópicas no están resueltas en detalle; tratarlas una por una exigiría un modelo de una escala completamente diferente. Paul Langevin, físico francés, propuso en 1908 representar su efecto agregado mediante una fuerza fluctuante junto a la fricción. El modelo resultante combina una evolución sistemática con una perturbación aleatoria cuya intensidad y estructura temporal deben especificarse.

La misma estrategia aparece en una población sometida a acontecimientos demográficos, en un circuito eléctrico con ruido térmico o en una variable fisiológica expuesta a perturbaciones continuas. En cada caso, sin embargo, debe justificarse qué representa la variable, de dónde provienen las fluctuaciones y cuándo resulta aceptable aproximarlas mediante ruido de pequeña escala temporal.

Si la fuerza restauradora es aproximadamente lineal en el desplazamiento respecto de un equilibrio, la relajación determinista tiene un tiempo característico $1/\kappa$ . Al incorporar incrementos brownianos, esa escala reaparece en la autocorrelación, mientras la intensidad del ruido determina el tamaño de las fluctuaciones. Los dos parámetros pueden estimarse a partir de propiedades distintas de una serie temporal, siempre que la observación sea compatible con la hipótesis de un proceso lineal estable.

La solución analítica es también un control para la integración numérica, porque en un paso de duración $h$ tanto la atenuación de la condición inicial como la varianza del ruido acumulado tienen expresiones exactas. Compararlas con Euler–Maruyama permite detectar errores que una trayectoria visualmente plausible podría ocultar.

131.2. Qué significa añadir ruido a una ecuación diferencial

La escritura informal

\frac{dX}{dt}=-\kappa(X-\mu)+\xi(t)

parece una ecuación diferencial ordinaria con un término adicional. Sin definir $\xi$ , sin embargo, no especifica un proceso matemático. Podría representar una fuerza periódica desconocida, un conjunto de impulsos discretos o una señal aleatoria correlacionada. También podría ser el llamado ruido blanco, que no existe como función ordinaria bien definida en cada instante. Para trabajar con este último caso necesitamos cambiar la notación.

Definiremos $W_t$ como un movimiento browniano estándar o proceso de Wiener. Comienza en cero, tiene trayectorias continuas, incrementos independientes en intervalos disjuntos y satisface

W_{t+h}-W_t\sim\mathcal N(0,h),\qquad h>0.

Por tanto, el incremento browniano típico es de orden $\sqrt h$ , mientras el incremento generado por una velocidad finita es de orden $h$ . Esta diferencia impide tratar $dW_t$ como una diferencial ordinaria. Una ecuación estocástica con ruido aditivo puede escribirse

\boxed{dX_t=-\kappa(X_t-\mu)\,dt+\sqrt{2D}\,dW_t},

donde $\kappa>0$ es la tasa de recuperación, $\mu$ el nivel de equilibrio y $D>0$ la intensidad del ruido. Llamaremos a este modelo proceso de Ornstein–Uhlenbeck, por los físicos George Uhlenbeck y Leonard Ornstein, que lo estudiaron en 1930 al analizar el movimiento browniano.

La ecuación dice que, durante un intervalo breve, hay una contribución sistemática aproximadamente igual a $-\kappa(X_t-\mu)h$ y una contribución aleatoria de desviación estándar $\sqrt{2Dh}$ . No predice el incremento concreto antes de observarlo. Define una distribución condicional para la evolución y, junto con la ley del estado inicial, determina la distribución de las trayectorias completas.

Las unidades ayudan a detectar errores. Si $X$ se mide en unidades $U$ y el tiempo en $T$ , entonces $\kappa$ tiene unidades $T^{-1}$ , $D$ tiene unidades $U^2/T$ y $\sqrt{2D}\,dW$ tiene unidades $U$ . A menudo se denomina «intensidad de ruido» tanto a $D$ como a $\sqrt{2D}$ , pese a que tienen unidades diferentes. Conviene explicitar la parametrización y las unidades, sobre todo al comparar artículos o programas.

Una trayectoria browniana es continua, pero casi seguramente no diferenciable en el sentido clásico. De ahí que su «velocidad instantánea de ruido» no pueda obtenerse mediante una derivada ordinaria. La expresión formal $\xi(t)=\sqrt{2D}\,dW_t/dt$ es útil como abreviatura de un ruido blanco generalizado, siempre que se recuerde esta limitación. En particular, la varianza de un ruido blanco ideal evaluado en un instante no es una magnitud finita comparable con la de una variable aleatoria ordinaria.

En el caso de ruido aditivo que estudiaremos, el coeficiente de $dW_t$ es constante, de modo que las interpretaciones de Itô y Stratonovich dan aquí la misma solución. Cuando la intensidad depende del propio estado, la elección de interpretación puede modificar la deriva efectiva. C133 y C134 se dedicarán a esa cuestión; por ahora bastará con entender qué modelo de incrementos estamos definiendo.

Once trayectorias de Ornstein–Uhlenbeck fluctúan en torno a una curva media que decae hacia 10 unidades. Una banda sombreada representa una desviación estándar teórica a cada tiempo. — Figura 131.1. Once realizaciones exactas del proceso de Ornstein–Uhlenbeck desde X₀=12 U. La curva central y la banda de ±1 desviación estándar representan estadísticas condicionales de conjunto, no un intervalo de confianza de la media. Semilla 13101. Ampliar figura ↗

131.3. Solución, dispersión y autocorrelación

La linealidad del proceso permite obtener una solución explícita. Escribamos $Y_t=X_t-\mu$ , de modo que $dY_t=-\kappa Y_t\,dt+\sqrt{2D}\,dW_t$ . Multiplicando por el factor integrante $e^{\kappa t}$ e integrando, se obtiene

Y_t=Y_0e^{-\kappa t} +\sqrt{2D}\int_0^t e^{-\kappa(t-s)}\,dW_s.

La primera parte es la relajación del desplazamiento inicial. La segunda reúne fluctuaciones ocurridas en momentos anteriores, con un peso que disminuye exponencialmente cuanto más lejano está el instante de la perturbación. Esta integral es estocástica; su integrador son los incrementos de Wiener.

Si comenzamos en un valor fijo $X_0=x_0$ , el término integral tiene esperanza cero. Por ello,

\boxed{E[X_t\mid X_0=x_0]=\mu+(x_0-\mu)e^{-\kappa t}}.

La media sigue exactamente la trayectoria del sistema determinista correspondiente. Las trayectorias individuales, en cambio, pueden cruzar repetidamente el equilibrio, alejarse de él o permanecer un tiempo en un mismo lado; solo la posición media del conjunto de realizaciones preparadas con el mismo estado inicial regresa de forma exponencial.

Para determinar su dispersión usamos la isometría de Itô, que permite calcular la varianza de una integral estocástica con integrando determinista. En este ejemplo,

\begin{aligned} \operatorname{Var}(X_t\mid X_0=x_0) &=2D\int_0^t e^{-2\kappa(t-s)}\,ds\\ &=\frac{D}{\kappa}\bigl(1-e^{-2\kappa t}\bigr). \end{aligned}

La varianza parte de cero porque hemos fijado el estado inicial exactamente. Crece de forma gradual y tiende a $D/\kappa$ cuando $t$ aumenta. En ausencia de restauración, $\kappa=0$ , tendríamos $X_t=X_0+\sqrt{2D}W_t$ y la varianza crecería como $2Dt$ ; para $\kappa>0$ , el retorno impide ese crecimiento indefinido.

Debido a que la integral de un integrando determinista frente a Wiener es gaussiana, la distribución condicional completa es

X_t\mid X_0=x_0\sim \mathcal N\!\left( \mu+(x_0-\mu)e^{-\kappa t}, \frac{D}{\kappa}(1-e^{-2\kappa t}) \right).

En esta expresión, el segundo argumento de la distribución normal es la varianza. Con la distribución condicional se calculan probabilidades de excedencia de umbrales o de intervalos, que la ecuación para la media no da.

Cuando $t$ crece, obtenemos una distribución límite normal con media $\mu$ y varianza $D/\kappa$ . Si escogemos $X_0$ aleatoriamente según esa distribución, de manera independiente de los incrementos brownianos futuros, el proceso es estacionario desde el comienzo. Si partimos de un $x_0$ fijo, inicialmente no lo es, aunque su distribución converja hacia la misma ley. Estabilidad determinista, convergencia en distribución y estacionariedad estricta son propiedades relacionadas, pero requieren enunciados distintos.

Los parámetros $\kappa$ y $D$ intervienen de maneras diferentes. Aumentar $D$ ensancha la distribución sin alterar la tasa de recuperación de la media. Reducir $\kappa$ aumenta simultáneamente la dispersión estacionaria y el tiempo de recuperación. Por eso una varianza amplia es compatible con perturbaciones fuertes, con recuperación lenta o con una combinación de ambas; distinguirlas exige información sobre la dinámica temporal.

Supongamos que el proceso ha sido inicializado en su distribución estacionaria. Su esperanza permanece en $\mu$ y su varianza es constante, aunque las observaciones sucesivas no son independientes. Un desplazamiento positivo suele continuar siendo positivo poco tiempo después, porque la restauración es gradual. Para cuantificar esa dependencia calcularemos la función de autocovarianza. Si $\tau$ es una separación temporal,

C_X(\tau) =E[(X_t-\mu)(X_{t+\tau}-\mu)].

La solución condicional permite escribir, para $\tau\geq0$ , $X_{t+\tau}-\mu=e^{-\kappa\tau}(X_t-\mu)+Z$ , donde $Z$ reúne perturbaciones posteriores a $t$ y es independiente de $X_t$ . El término cruzado tiene esperanza cero. En consecuencia,

\boxed{C_X(\tau)=\frac{D}{\kappa}e^{-\kappa|\tau|}}, \qquad \rho_X(\tau)=e^{-\kappa|\tau|}.

La autocorrelación normalizada $\rho_X$ permite distinguir la persistencia temporal de la amplitud de las fluctuaciones. El parámetro $D$ determina la escala vertical de la autocovarianza; $\kappa$ , la velocidad con que se pierde la correlación. Definimos el tiempo de correlación integral, para este proceso de correlación positiva, como

\tau_{\mathrm{corr}}=\int_0^{\infty}\rho_X(s)\,ds=\frac1\kappa.

La misma escala aparece en la relajación determinista. Esa coincidencia es propia de esta ecuación lineal. En un sistema con diferentes modos de relajación podrían coexistir varias escalas temporales y las correlaciones podrían cambiar de signo. Una función de autocorrelación observada también puede resultar de la combinación de diversos mecanismos.

El espectro de potencia es otra representación de la misma dependencia. Usaremos una densidad espectral bilateral con frecuencia angular $\omega$ , definida mediante la transformada de Fourier

S_X(\omega)=\int_{-\infty}^{\infty}C_X(\tau)e^{-i\omega\tau}\,d\tau.

Para la autocovarianza exponencial se obtiene

\boxed{S_X(\omega)=\frac{2D}{\kappa^2+\omega^2}}.

El espectro es más intenso a frecuencias bajas y disminuye cuando $|\omega|$ aumenta. Aquí la variable observada $X$ es coloreada, aunque el término que la impulsa sea ruido blanco. Así, a partir de la correlación de una salida no puede deducirse directamente la del estímulo que la genera, porque el sistema transforma las propiedades espectrales de sus entradas.

Hay además una consecuencia para la inferencia estadística. Si registramos mil puntos de una trayectoria con intervalos mucho más cortos que $1/\kappa$ , tenemos mil mediciones, pero los puntos cercanos contienen información parcialmente redundante. Una estimación de la media que tratara todos los puntos como independientes daría, en general, una precisión excesivamente optimista. La precisión inferencial depende de la autocorrelación de la serie, de la duración del registro y del modelo usado para el ruido.

Seis densidades normales: cinco distribuciones transitorias desde el estado inicial fijo 12 y la distribución estacionaria de media 10. La posición media converge y la dispersión aumenta. — Figura 131.2. Distribuciones condicionales exactas del proceso restaurador en cinco tiempos y su límite estacionario N(10,0,5). Las curvas representan leyes distintas en cada instante; no son historiales de una trayectoria. Ampliar figura ↗

131.4. Movimiento browniano, fricción y fluctuación–disipación

Una de las primeras aplicaciones del proceso de Ornstein–Uhlenbeck fue la velocidad de una partícula browniana. Sea $m$ su masa, $v(t)$ una componente de velocidad y $\gamma$ el coeficiente de fricción lineal, definido mediante la fuerza $-\gamma v$ . La ecuación de Langevin idealizada correspondiente se escribe

m\,dv_t=-\gamma v_t\,dt+\sqrt{2\gamma k_B T}\,dW_t,

donde $k_B$ es la constante de Boltzmann y $T$ la temperatura absoluta del baño térmico. Se supone un entorno isotermo en equilibrio, fricción lineal y ruido térmico idealizado sin correlación temporal resuelta. Estas hipótesis restringen el modelo; en fluidos con memoria hidrodinámica, por ejemplo, la fricción instantánea puede requerir un modelo generalizado con un núcleo temporal.

Dividiendo por la masa identificamos una tasa de relajación $\kappa=\gamma/m$ y una intensidad de difusión en el espacio de velocidades

D_v=\frac{\gamma k_B T}{m^2}.

La varianza estacionaria de una componente de velocidad es

\operatorname{Var}(v)=\frac{D_v}{\kappa}=\frac{k_B T}{m}.

Así recuperamos $mE[v^2]/2=k_B T/2$ , la equipartición clásica de energía para ese grado de libertad. El factor que relaciona la fricción con la amplitud del ruido viene impuesto por la consistencia con el equilibrio térmico. En este ejemplo aparece una forma de la relación entre fluctuación y disipación, que Ryogo Kubo estudió después en términos generales.

Podemos pasar de velocidad a posición. En la aproximación sobreamortiguada, cuando examinamos tiempos bastante mayores que $m/\gamma$ y la inercia puede despreciarse para el fenómeno de interés, una partícula en un potencial $U(x)$ obedece

dx_t=-\frac{U'(x_t)}{\gamma}\,dt +\sqrt{\frac{2k_BT}{\gamma}}\,dW_t.

Si el potencial es armónico, $U(x)=kx^2/2$ , la ecuación vuelve a ser de Ornstein–Uhlenbeck con $\kappa=k/\gamma$ y $D_x=k_BT/\gamma$ . Su distribución estacionaria tiene

\operatorname{Var}(x)=\frac{k_BT}{k},\qquad \tau_{\mathrm{corr}}=\frac{\gamma}{k}.

Una trampa más rígida reduce las fluctuaciones de posición, mientras una fricción mayor enlentece la recuperación. En este caso particular la amplitud estacionaria no depende de $\gamma$ , porque la disipación y la excitación térmica cambian conjuntamente. Si variáramos la fricción manteniendo artificialmente fija una intensidad de ruido externa, ese resultado ya no se conservaría.

Esta derivación también marca un límite. El teorema de fluctuación–disipación descansa en supuestos físicos de equilibrio y de respuesta lineal, que las series económicas, poblacionales o clínicas no tienen por qué cumplir. En esos sistemas puede usarse una ecuación parecida para describir fluctuaciones, pero la razón entre disipación e intensidad de ruido debe estimarse o justificarse por separado. La semejanza formal de las ecuaciones no es evidencia de una temperatura efectiva con sentido físico.

A la izquierda, densidades gaussianas de velocidad para tres temperaturas. A la derecha, la varianza teórica igual a la temperatura reducida y valores empíricos de simulaciones exactas. — Figura 131.3. En el modelo lineal de Langevin para la velocidad, la relación entre fricción y ruido térmico produce var(v)=kB T/m. Las 30 000 realizaciones por temperatura se obtuvieron con la transición exacta a t=10 desde v₀=0; son datos simulados, no mediciones de partículas. Ampliar figura ↗

131.5. Las ecuaciones de momentos y la relación con una distribución

Los resultados anteriores se obtuvieron resolviendo la ecuación completa. Conviene tener también una vía que no requiera conocer de antemano la densidad de probabilidad. La ecuación diferencial estocástica establece cómo se mueve cada realización; la evolución de sus momentos resume ciertas propiedades de la distribución. Para el proceso lineal, el primer momento satisface

\frac{d}{dt}E[X_t]=-\kappa(E[X_t]-\mu).

Esto se obtiene tomando esperanzas en la ecuación integral, porque el incremento estocástico tiene esperanza cero. La ecuación se cierra en la propia media, sin que intervengan la varianza ni los momentos superiores. Ese cierre depende de que la deriva sea lineal. En modelos con deriva no lineal, $E[f(X)]$ normalmente no coincide con $f(E[X])$ y aparecen jerarquías de ecuaciones acopladas.

El segundo momento requiere tener en cuenta una propiedad específica de la fluctuación browniana. Escribiendo $Y=X-\mu$ , el cálculo estocástico da

d(Y_t^2)=2Y_t\,dY_t+(dY_t)^2,

con la regla $(dW_t)^2=dt$ en el sentido de variación cuadrática. Como el coeficiente del ruido es $\sqrt{2D}$ , el término adicional vale $2D\,dt$ . Al tomar esperanzas queda

\frac{d}{dt}E[Y_t^2]=-2\kappa E[Y_t^2]+2D.

Esta relación puede verificarse directamente desde la solución exacta; aquí se muestra para anticipar por qué una regla de cálculo ordinaria omitiría una contribución esencial. El capítulo C133 establecerá formalmente la regla de Itô y su campo de aplicación.

Si $m_Y(t)=E[Y_t]$ , la varianza $V(t)=E[Y_t^2]-m_Y(t)^2$ satisface

\boxed{\frac{dV}{dt}=-2\kappa V+2D}.

Su solución general es

V(t)=V(0)e^{-2\kappa t}+\frac D\kappa(1-e^{-2\kappa t}).

Esta expresión admite una condición inicial aleatoria de varianza finita; el resultado del apartado 131.3 correspondía al caso especial $V(0)=0$ . La varianza estacionaria es el equilibrio entre una contribución que la reduce, proporcional a $2\kappa V$ , y otra que la alimenta, de magnitud $2D$ . Esta interpretación es válida en el nivel de los momentos y no convierte el ruido blanco en una serie de «empujones» observables individualmente.

También podemos describir la evolución de la densidad $p(x,t)$ . Bajo condiciones regulares, la ecuación que gobierna esa densidad es

\frac{\partial p}{\partial t}= \kappa\frac{\partial}{\partial x}[(x-\mu)p] +D\frac{\partial^2p}{\partial x^2}.

Esta es la ecuación de Fokker–Planck para el proceso considerado. Su lado derecho combina transporte probabilístico generado por la deriva y difusión en el espacio de estados. Aquí basta con constatar que el mismo modelo admite una descripción en términos de trayectorias y otra en términos de densidades, cada una adecuada para preguntas diferentes. En C135 desarrollaremos la ecuación de Fokker–Planck con más detalle, incluyendo corrientes de probabilidad y condiciones de frontera. La densidad no recorre una trayectoria individual; describe la ley de un conjunto de realizaciones.

131.6. Recuperación lenta, fluctuaciones grandes y señales de alerta temprana

Un resultado del modelo tiene interés fuera del movimiento browniano. La varianza estacionaria es $D/\kappa$ y el tiempo de correlación es $1/\kappa$ . Si la recuperación del equilibrio se hace más lenta, $\kappa$ disminuye; con ruido aproximadamente constante, aumentan tanto la dispersión como la persistencia de las fluctuaciones. Esta relación aparece con frecuencia en discusiones sobre critical slowing down, o enlentecimiento crítico, cerca de ciertas transiciones dinámicas.

La intuición puede examinarse mediante un modelo no lineal $\dot x=f(x;r)$ , con parámetro de control $r$ y equilibrio $x^*(r)$ . Si $f(x^*;r)=0$ , cerca del equilibrio una pequeña perturbación $y=x-x^*$ satisface aproximadamente

\dot y\approx f_x(x^*;r)y.

Cuando $f_x(x^*;r)<0$ , identificamos $\kappa=-f_x(x^*;r)$ . Si un autovalor real de estabilidad se aproxima a cero por valores negativos, la recuperación lineal se hace lenta. Bajo un forzamiento estocástico adecuado, algunas estadísticas pueden reflejar ese cambio. La interpretación requiere que el sistema esté suficientemente cerca del equilibrio para justificar la linealización y que las fluctuaciones sean pequeñas en comparación con la escala donde intervienen los términos no lineales.

No todas las bifurcaciones comparten la misma firma. En un sistema multidimensional pueden intervenir modos oscilatorios, correlaciones cruzadas y proyecciones observacionales que ocultan el modo lento. Una variable medida podría mostrar cambios modestos incluso si otra concentra la pérdida de estabilidad. La velocidad con que varía el parámetro de control determina también si resulta apropiado aproximar cada ventana temporal mediante una distribución estacionaria. Si el sistema atraviesa rápidamente una transición, las estadísticas observadas pueden reflejar retrasos y transitorios antes que propiedades de equilibrio local.

En una serie real, el aumento de varianza o de autocorrelación también puede tener otras causas. Quizá aumentó la intensidad de las perturbaciones, cambió la resolución del instrumento, aparecieron tendencias externas o se modificó la composición de la población muestreada. La correlación temporal estimada en ventanas breves tiene incertidumbre considerable. Una tendencia aparente puede incluso surgir al seleccionar retrospectivamente las ventanas anteriores a un acontecimiento.

Las llamadas señales de alerta temprana son, en consecuencia, hipótesis contrastables sobre la proximidad de una transición. El modelo da un mecanismo que predice relaciones específicas entre varianza, tiempo de correlación y tasa local de recuperación; para atribuir esas relaciones a una transición concreta hacen falta hipótesis rivales, estimaciones de incertidumbre y conocimiento del proceso observado.

Dos curvas muestran que al disminuir la tasa restauradora κ aumenta la varianza D/κ y crece la autocorrelación al mismo retardo, manteniendo D fija. — Figura 131.4. Predicciones exactas del modelo lineal estacionario con D fija: var(X)=D/κ y ρ(Δ)=e^(−κΔ), para Δ=1 min. La asociación no constituye por sí sola una señal validada de bifurcación ni implica que los datos reales sean estacionarios. Ampliar figura ↗

131.7. Simulación e inferencia a partir de series

La forma habitual de aproximar la ecuación de Langevin en intervalos $\Delta t$ es el método de Euler–Maruyama. Para el modelo lineal,

X_{n+1}=X_n-\kappa(X_n-\mu)\Delta t +\sqrt{2D\Delta t}\,Z_n, \qquad Z_n\sim\mathcal N(0,1),

con $Z_n$ independientes. Esta regla aproxima la dinámica continua, pero también define una cadena discreta distinta. El factor de recuperación de su media es $1-\kappa\Delta t$ ; el del proceso continuo, $e^{-\kappa\Delta t}$ . Si el paso es demasiado grande puede introducir oscilaciones artificiales o amplificar los errores. La estabilidad de su segundo momento exige $|1-\kappa\Delta t|<1$ , equivalente a $0<\kappa\Delta t<2$ .

El sesgo puede calcularse sin ejecutar un programa. Si la condición anterior se cumple, la varianza estacionaria de la discretización es

V_{\mathrm{EM}}=\frac{2D\Delta t}{1-(1-\kappa\Delta t)^2} =\frac{2D}{2\kappa-\kappa^2\Delta t}.

Para $\Delta t>0$ , esta expresión supera $D/\kappa$ . La simulación puede mostrar fluctuaciones excesivas aunque el algoritmo parezca estable. En un modelo no lineal, cuyos momentos no admitan expresiones cerradas, el efecto podría ser más difícil de descubrir.

El proceso de Ornstein–Uhlenbeck tiene la ventaja de que puede simularse en tiempos igualmente espaciados sin error de discretización temporal. La distribución de transición exacta es

\boxed{X_{n+1}=\mu+e^{-\kappa\Delta t}(X_n-\mu) +\sqrt{\frac D\kappa(1-e^{-2\kappa\Delta t})}\,Z_n}.

Aquí el paso determina solamente los instantes a los que observamos la trayectoria; la distribución conjunta de esos valores coincide con la del proceso continuo. Los movimientos entre muestras siguen siendo desconocidos y no deben interpretarse como segmentos rectos reales. Para un modelo con coeficientes no lineales normalmente no contaremos con esta actualización exacta, pero el caso lineal es una referencia para verificar métodos aproximados.

El siguiente programa compara las distribuciones estacionarias producidas por ambos procedimientos. Usa muchas trayectorias independientes y un horizonte suficientemente largo respecto de $1/\kappa$ . La semilla permite repetir el experimento, mientras los controles numéricos detectan discrepancias importantes con los valores analíticos.

import numpy as np

rng = np.random.default_rng(131)
kappa, D, mu = 1.2, 0.45, 0.4
dt, n_steps, n_paths = 0.2, 80, 60_000
x_exact = np.full(n_paths, 2.4)
x_euler = x_exact.copy()
a = np.exp(-kappa * dt)
s_exact = np.sqrt((D / kappa) * (1 - a*a))
s_euler = np.sqrt(2 * D * dt)

for _ in range(n_steps):
    z = rng.standard_normal(n_paths)
    x_exact = mu + a * (x_exact - mu) + s_exact * z
    x_euler = x_euler - kappa * (x_euler - mu) * dt + s_euler * z

var_real = D / kappa
var_euler = 2 * D / (2*kappa - kappa*kappa*dt)
print('Media exacta estimada:', round(x_exact.mean(), 4))
print('Varianza exacta estimada / teórica:',
      round(x_exact.var(), 4), round(var_real, 4))
print('Varianza Euler estimada / teórica:',
      round(x_euler.var(), 4), round(var_euler, 4))

assert abs(x_exact.mean() - mu) < 0.02
assert abs(x_exact.var() - var_real) < 0.02
assert abs(x_euler.var() - var_euler) < 0.025
assert var_euler > var_real

La actualización exacta aproxima la distribución estacionaria del proceso continuo a partir de valores observados en tiempos discretos. Euler–Maruyama converge hacia ella cuando se reduce el paso, pero conserva un sesgo para $\Delta t$ finito. Las tolerancias indicadas son controles de funcionamiento y no reemplazan un análisis de error Monte Carlo; pueden estrecharse aumentando el número de realizaciones y usando réplicas independientes. Si una simulación de un modelo desconocido produce un fenómeno sorprendente, hay que comprobar que persiste al variar la resolución temporal.

Otra forma de auditar la simulación es comparar las transiciones condicionales antes de alcanzar la distribución estacionaria. Si todas las realizaciones empiezan en un mismo punto, la teoría da la media decreciente y la varianza creciente en cada instante. Una prueba que verificara solamente el resultado final podría pasar por alto un algoritmo que alcanza aproximadamente la distribución correcta mediante una trayectoria transitoria equivocada. Conviene medir ambas magnitudes en varios tiempos, en particular cerca de una constante de relajación y durante el acercamiento al régimen estacionario.

También hay que distinguir entre exactitud de las marginales y exactitud del proceso completo. Un generador que produjera muestras gaussianas independientes con media $\mu$ y varianza $D/\kappa$ reproduciría perfectamente el histograma estacionario. Sin embargo, destruiría la correlación temporal $e^{-\kappa|\tau|}$ y daría resultados incorrectos sobre episodios persistentes, velocidades de recuperación inferidas y cruces consecutivos de umbrales. Validar un proceso estocástico exige entonces comprobar también las dependencias entre varios tiempos cuando forman parte de las predicciones del modelo.

Finalmente, una condición de frontera puede modificar drásticamente la interpretación. El proceso de Ornstein–Uhlenbeck sin barreras puede tomar cualquier valor real, aunque los extremos sean poco probables. Si el fenómeno contiene una barrera absorbente, un techo físico o una regla de intervención al cruzar cierto umbral, la evolución posterior ya no coincide con la solución libre. El tiempo necesario para alcanzar por primera vez el umbral tampoco se deduce de una sola probabilidad marginal de excedencia. Los problemas de primera llegada y escape se estudiarán con herramientas específicas en C137 y C138. Introducir límites físicos como simples recortes al dibujar un gráfico escondería una modificación sustantiva del modelo.

Imaginemos que registramos una variable de un sistema aparentemente estable cada cinco minutos. La señal oscila alrededor de una referencia y presenta cierta autocorrelación. Un proceso de Ornstein–Uhlenbeck podría ser una primera descripción, pero ajustarlo requiere distinguir la variación del estado, las propiedades de la perturbación y el error de medición. La formulación precisa dice que, condicionado a un estado conocido $X_t=x$ , la media y la varianza del estado transcurrido un intervalo $\Delta$ son

E[X_{t+\Delta}\mid X_t=x] =\mu+e^{-\kappa\Delta}(x-\mu),

\operatorname{Var}(X_{t+\Delta}\mid X_t=x) =\frac D\kappa(1-e^{-2\kappa\Delta}).

Estas relaciones permiten construir una verosimilitud para observaciones discretas sin reemplazar la evolución exacta por Euler. Los residuos estandarizados del modelo deberían comportarse como innovaciones independientes si la observación es exacta y los supuestos se cumplen. Un patrón persistente de correlación residual sugiere que falta estructura temporal, que cambian los parámetros o que la ecuación de observación requiere ajustes.

La frecuencia de muestreo limita la identificación. Si $\Delta$ es mucho mayor que $1/\kappa$ , los estados observados estarán casi descorrelacionados; resultará difícil estimar con precisión la velocidad de recuperación porque las correlaciones informativas ocurren entre mediciones. Si $\Delta$ es muy pequeño en relación con esa escala y el registro dura poco, abundarán valores casi idénticos sin suficiente información sobre la dinámica lenta.

El ruido instrumental agrega otra dificultad. Si observamos $Y_n=X_{t_n}+\varepsilon_n$ , con errores de medición independientes, la varianza observada incluye una contribución extra. Las covarianzas entre tiempos distintos pueden conservar parte de la información del proceso latente, pero la autocorrelación normalizada de la señal medida cambia porque el denominador ahora incorpora error. Una aplicación ingenua podría atribuir al mecanismo físico una recuperación excesivamente rápida. Un modelo de espacio de estados permite distinguir, bajo supuestos identificables, evolución latente y medición.

También importa cómo se eligió el equilibrio. Restar una tendencia estimada a la propia serie puede eliminar fluctuaciones lentas genuinas; no restarla cuando existe deriva sistemática puede crear autocorrelaciones espurias. El equilibrio de un sistema físico y la media calculada durante una ventana son conceptos distintos. Llamar «estacionaria» a una serie exige evaluar si sus propiedades probabilísticas permanecen aproximadamente invariantes en el intervalo considerado.

En investigación biomédica conviene ser cauto al interpretar estos parámetros como rasgos clínicos. Una tasa de recuperación matemática describe el comportamiento de una variable bajo un protocolo, sin aludir a una capacidad psicológica general. Puede ser útil para comparar modelos o generar hipótesis, pero su interpretación causal exige identificar entradas, mecanismos y escalas. La aparente universalidad del proceso de Ornstein–Uhlenbeck proviene de la simplicidad de su ecuación lineal.

Curva creciente del sesgo positivo de la varianza estacionaria aproximada cuando aumenta κh. La integración exacta tiene sesgo cero. Puntos simulados reproducen el sesgo de Euler. — Figura 131.5. Para el proceso lineal dX=−κXdt+√(2D)dW, Euler–Maruyama produce var_EM=(D/κ)/(1−κh/2), si 0<κh<2. Los puntos son realizaciones de la recurrencia gaussiana de Euler después de 400 pasos; la transición exacta conserva D/κ para cualquier h>0. Las divergencias próximas al límite de estabilidad hacen inadecuados los pasos grandes. Ampliar figura ↗

131.8. Comprobación del modelo lineal

Consideremos una variable registrada en unidades arbitrarias, con equilibrio $\mu=10$ , tasa de recuperación $\kappa=0{,}4\ \mathrm{min}^{-1}$ e intensidad $D=0{,}2\ U^2/\mathrm{min}$ . Una observación inicial es $X_0=12$ . La distancia media respecto del equilibrio, al cabo de cinco minutos, será

E[X_5]-10=2e^{-0{,}4\cdot5}=2e^{-2}\approx0{,}271.

Por tanto, el valor esperado del conjunto de realizaciones que parten de ese valor inicial será aproximadamente $10{,}271$ . La varianza condicional a los cinco minutos será

\operatorname{Var}(X_5\mid X_0=12) =\frac{0{,}2}{0{,}4}(1-e^{-4}) \approx0{,}4908\ U^2.

La desviación estándar condicional es, por tanto, cercana a $0{,}701\ U$ . La distribución estacionaria tiene media 10, varianza $0{,}5\ U^2$ y desviación estándar aproximadamente $0{,}707\ U$ . El tiempo de correlación característico es $2{,}5$ minutos, de modo que dos observaciones estacionarias separadas cinco minutos tienen correlación $e^{-2}\approx0{,}1353$ , bajo los supuestos del modelo.

Estas cantidades responden a preguntas distintas. La media condicional cuantifica la recuperación respecto de un estado inicial; la varianza condicional describe la incertidumbre entre realizaciones al tiempo considerado; la autocorrelación estacionaria mide la persistencia temporal cuando el proceso se encuentra en equilibrio estadístico. Mezclarlas conduce a errores habituales, como calcular la precisión de una media experimental usando un número de observaciones independientes que en realidad no se tiene.

El proceso de Ornstein–Uhlenbeck permite separar el tiempo de recuperación de la dispersión sostenida por las fluctuaciones. La tasa $\kappa$ controla el decaimiento de la dependencia temporal y la razón $D/\kappa$ fija la varianza estacionaria. La comparación con Euler–Maruyama muestra que un paso temporal grande puede inflar esa varianza aun cuando el algoritmo conserve una media razonable.

En aplicaciones físicas, la relación entre amplitud de ruido y disipación requiere condiciones termodinámicas precisas; una ecuación restauradora cualquiera no cumple por ello una relación de fluctuación–disipación de equilibrio. La versión exacta de tiempo discreto es un contraste numérico especialmente útil cuando el objetivo es estimar correlaciones o analizar datos muestreados a intervalos fijos.

Referencias

Gardiner, C. W. (2009). Stochastic Methods: A Handbook for the Natural and Social Sciences (4.ª ed.). Springer. Procesos de difusión, Ornstein–Uhlenbeck y cálculo de correlaciones.
Kubo, R. (1966). The fluctuation-dissipation theorem. Reports on Progress in Physics, 29(1), 255–284. Relación entre respuesta, disipación y fluctuaciones en equilibrio. DOI: 10.1088/0034-4885/29/1/306.
Langevin, P. (1908). Sur la théorie du mouvement brownien. Comptes Rendus de l'Académie des Sciences, 146, 530–533. Formulación histórica del modelo de fricción y fuerza fluctuante.
Risken, H. (1989). The Fokker–Planck Equation: Methods of Solution and Applications (2.ª ed.). Springer. Dinámica de densidades y procesos lineales.
Uhlenbeck, G. E., y Ornstein, L. S. (1930). On the theory of the Brownian motion. Physical Review, 36(5), 823–841. Descripción de la dinámica de velocidad browniana. DOI: 10.1103/PhysRev.36.823.
van Kampen, N. G. (2007). Stochastic Processes in Physics and Chemistry (3.ª ed.). Elsevier. Aproximaciones de Langevin, ecuaciones maestras y relación entre niveles de descripción.