Euler–Maruyama y verificación de una simulación estocástica

Año 2 · Módulo 14

140.1. Verificación de trayectorias numéricas

Un integrador estocástico puede reproducir la media final de un proceso y, al mismo tiempo, aproximar mal sus trayectorias. Otro puede representar adecuadamente las trayectorias y fallar en la probabilidad de un evento raro. Estas diferencias obligan a definir el criterio de exactitud antes de elegir el paso numérico. Euler–Maruyama, que usa incrementos brownianos de varianza proporcional al intervalo temporal, es el punto de partida para ese análisis y permite estudiar convergencia fuerte, convergencia débil y sesgos de frontera.

En una ecuación diferencial ordinaria, un algoritmo aproxima una trayectoria determinada por condiciones iniciales y parámetros. En una ecuación diferencial estocástica también tenemos una trayectoria, pero definida conjuntamente por el estado inicial y una realización de la fuente aleatoria. Si comparamos dos ejecuciones con semillas distintas, las diferencias pueden deberse a trayectorias brownianas distintas, aun cuando ambos algoritmos sean correctos. Para evaluar el error temporal necesitamos distinguir la aproximación de una misma realización y la aproximación de propiedades estadísticas obtenidas al promediar realizaciones.

Consideremos una difusión de Itô en una dimensión

dX_t=a(X_t,t)\,dt+b(X_t,t)\,dW_t,\qquad X_0=x_0.

La deriva $a$ es la contribución sistemática al cambio y $b$ determina la amplitud del ruido. $W_t$ es un movimiento browniano estándar, con incrementos independientes de distribución normal y varianza igual a la duración del intervalo. Para interpretar el modelo deben conocerse el dominio de $X$ , las unidades de los coeficientes y la condición inicial. También importa decidir si las fronteras son reflectantes, absorbentes, periódicas o simplemente límites numéricos sin significado físico.

El método de Euler–Maruyama construye una sucesión discreta que aproxima $X_t$ en instantes seleccionados. Sus autores epónimos son Leonhard Euler, matemático del siglo XVIII, y Gisiro Maruyama, matemático japonés que trabajó sobre ecuaciones diferenciales estocásticas en el siglo XX. El algoritmo comparte con Euler la aproximación de primer paso para la deriva, pero introduce incrementos brownianos para representar la integral estocástica. Es una herramienta de referencia por su simplicidad y porque permite establecer un protocolo reproducible de verificación.

La dificultad principal consiste en no exigirle lo que no garantiza. En un problema, quizás nos interese la distribución de $X_T$ ; en otro, el tiempo en que $X_t$ cruza un umbral; en otro, la probabilidad de supervivencia de un sistema. Esas cantidades responden de manera diferente a la discretización. Esta clase desarrolla las actualizaciones de Euler–Maruyama, sus órdenes de convergencia y sus fallos conocidos, junto con pruebas que permiten detectar errores antes de publicar conclusiones.

140.2. La malla temporal y el tamaño correcto del incremento aleatorio

Definamos una partición regular $t_n=nh$ , con $n=0,\ldots,N$ , $Nh=T$ y $h>0$ . El paso $h$ tiene unidades de tiempo. La propiedad fundamental del movimiento browniano es

\Delta W_n=W_{t_{n+1}}-W_{t_n}\sim\mathcal N(0,h).

Podemos generar cada incremento como $\Delta W_n=\sqrt h\,Z_n$ , donde los $Z_n$ son independientes y siguen una normal estándar. Aparece una raíz cuadrada donde las integrales ordinarias tendrían el factor $h$ . El tamaño típico del movimiento aleatorio en un intervalo pequeño escala como $\sqrt h$ , aunque su valor esperado sea cero.

Un error elemental de graves consecuencias es usar $hZ_n$ en vez de $\sqrt hZ_n$ . Si hiciéramos eso, la varianza acumulada durante $T$ sería proporcional a $Nh^2=Th$ y tendería a cero al refinar la malla. Al mejorar la resolución estaríamos eliminando artificialmente las fluctuaciones que la ecuación especifica.

La forma integral de la ecuación es

X_{t_{n+1}}=X_{t_n} +\int_{t_n}^{t_{n+1}}a(X_s,s)\,ds +\int_{t_n}^{t_{n+1}}b(X_s,s)\,dW_s.

El método aproxima los coeficientes durante el intervalo mediante sus valores en el extremo izquierdo. La actualización queda

\boxed{X_{n+1}=X_n+a(X_n,t_n)h+b(X_n,t_n)\sqrt h\,Z_n.}

Se trata de una discretización de Itô. Evaluar el coeficiente de ruido en el extremo izquierdo corresponde a la convención de integración que adopta el modelo. Si sustituyéramos esa evaluación por el punto medio sin compensar los términos de deriva, podríamos aproximar otro proceso. La clase C134 mostró las condiciones bajo las cuales se convierten correctamente los coeficientes entre Itô y Stratonovich.

La independencia de los incrementos debe respetarse entre intervalos no superpuestos de una misma trayectoria. Si simulamos múltiples realizaciones independientes, también necesitaremos fuentes aleatorias independientes entre ellas, salvo que usemos deliberadamente un acoplamiento para comparar algoritmos. La semilla de un generador se conserva para reproducibilidad, pero no es una característica física de la ecuación. Cambiarla modifica los valores muestrales; una conclusión robusta debe resistir ese cambio.

La malla temporal también tiene un coste. Reducir $h$ a la mitad duplica el número de pasos para un horizonte fijo. Si el objetivo requiere millones de realizaciones o tiempos de integración largos, no es razonable prometer precisión arbitraria. La verificación establece qué error puede tolerarse frente a los recursos disponibles y qué magnitudes merecen un refinamiento específico.

Panel izquierdo con una trayectoria browniana y sus muestras en tres resoluciones; panel derecho con desviación típica proporcional a la raíz del paso. — Figura 140.1. Una realización de W(t) muestreada en mallas anidadas (512, 32 y 16 pasos); cada incremento grueso suma los incrementos finos. En el panel derecho la desviación típica de ΔW confirma la escala √h, frente a la actualización incorrecta de orden h. Semillas 14001 y 14011. Ampliar figura ↗

140.3. Derivación de Euler–Maruyama y papel de la variación cuadrática

¿Por qué no basta con aplicar el método de Euler ordinario y tratar el ruido como una función cualquiera del tiempo? La dificultad es que las trayectorias brownianas son continuas, pero casi seguramente no diferenciables en ningún punto. La expresión informal $dW_t/dt$ no es una derivada ordinaria. Su papel matemático se establece mediante una integral estocástica y sus propiedades de convergencia.

Para una deriva suficientemente regular, la integral temporal durante un intervalo corto satisface aproximadamente $\int_{t_n}^{t_{n+1}}a(X_s,s)ds\approx a(X_n,t_n)h$ . Para la parte estocástica, si el coeficiente puede considerarse constante en el mismo intervalo, obtenemos $\int_{t_n}^{t_{n+1}}b(X_s,s)dW_s\approx b(X_n,t_n)\Delta W_n$ . La suma de ambas contribuciones produce Euler–Maruyama. La precisión depende de cómo varían los coeficientes y el estado dentro del intervalo.

La diferencia respecto del cálculo ordinario aparece al examinar los términos omitidos. Para un incremento browniano, $(\Delta W_n)^2$ tiene esperanza $h$ , del mismo orden que el paso. En una partición cada vez más fina,

\sum_{n=0}^{N-1}(\Delta W_n)^2\longrightarrow T

en probabilidad; de hecho, converge también en otros sentidos bajo construcciones habituales. Esta variación cuadrática explica el término adicional de la fórmula de Itô. También explica por qué los órdenes de convergencia de los métodos estocásticos no se obtienen copiando directamente los del cálculo determinista.

Bajo hipótesis apropiadas, la expansión estocástica de orden superior conduce, para una difusión escalar, al método de Milstein

X_{n+1}=X_n+a(X_n)h+b(X_n)\Delta W_n +\frac12b(X_n)b'(X_n)\left[(\Delta W_n)^2-h\right].

El último término surge al conservar una contribución de la expansión de Itô que Euler–Maruyama descarta. Para un coeficiente de difusión constante, $b'(x)=0$ y esa corrección desaparece. En problemas con ruido multiplicativo, como el movimiento browniano geométrico, puede mejorar el orden fuerte de convergencia.

Esta expresión de Milstein presupone ruido escalar y coeficientes suficientemente regulares. La generalización a sistemas multidimensionales con varios ruidos puede exigir términos iterados adicionales. Sería incorrecto aplicar mecánicamente una fórmula unidimensional a un sistema en el que los campos de difusión no conmutan. Tampoco debemos suponer que la mera presencia del término corrector resuelve problemas de estabilidad, positividad o cruces de fronteras.

Desde el punto de vista computacional, la derivación justifica dos decisiones prácticas. La primera es registrar qué interpretación estocástica se usa, porque las mismas letras $a$ y $b$ pueden describir procesos diferentes según esa elección. La segunda es conservar los incrementos brownianos en los experimentos de convergencia, porque solo así se comparan dos discretizaciones de una misma trayectoria subyacente.

140.4. Errores fuertes y débiles como preguntas científicas diferentes

Supongamos que disponemos de una solución exacta $X_T$ y de una aproximación $X_N^{(h)}$ obtenida con la misma trayectoria browniana. El error fuerte mide la proximidad entre realizaciones acopladas. Una medida habitual es la raíz del error cuadrático medio

e_{\mathrm{fuerte}}(h)=\left(\mathbb E\left[|X_T-X_N^{(h)}|^2\right]\right)^{1/2}.

Si $e_{\mathrm{fuerte}}(h)\leq C h^p$ para pasos suficientemente pequeños, hablamos de orden fuerte $p$ en esa norma. Bajo condiciones estándar de Lipschitz global, crecimiento lineal y regularidad apropiada, Euler–Maruyama tiene típicamente orden fuerte $1/2$ . Esto significa que reducir el paso a la cuarta parte reduce aproximadamente a la mitad el error fuerte en el régimen asintótico; no garantiza ese factor para cualquier paso inicial ni cualquier sistema.

El error débil estudia expectativas de observables. Para una función de prueba $\varphi$ , consideramos

e_{\mathrm{débil}}(h;\varphi) =\left|\mathbb E[\varphi(X_T)]-\mathbb E[\varphi(X_N^{(h)})]\right|.

Con coeficientes y observables suficientemente suaves, Euler–Maruyama suele alcanzar orden débil uno. Una simulación puede, por tanto, estimar bastante bien una media mientras difiere apreciablemente de la realización exacta correspondiente. No hay contradicción, porque medimos propiedades distintas. El orden débil tampoco está garantizado para indicadores discontinuos, probabilidades de cruce o funciones singulares cerca de una frontera.

La distinción afecta al diseño de la prueba. Para estimar el error fuerte, se deben construir discretizaciones de distintas resoluciones sobre un mismo ruido subyacente. Una manera sencilla es generar todos los incrementos de una malla fina y sumar grupos consecutivos para formar los incrementos de las mallas gruesas. Si solo cambiamos $h$ y dejamos que cada ejecución produzca incrementos independientes, la diferencia entre resultados contiene una variabilidad que no desaparece como un error de aproximación de trayectoria.

Para el error débil interesa comparar medias u otras expectativas, pero una media obtenida a partir de $M$ trayectorias contiene error de Monte Carlo de orden $M^{-1/2}$ bajo supuestos habituales de varianza finita. Refinar $h$ no reduce esa incertidumbre estadística. Si el sesgo temporal es menor que el error muestral, la gráfica del error débil puede quedar dominada por el ruido de estimación y parecer errática.

Hay además una distinción entre error respecto del modelo matemático y adecuación del modelo al fenómeno. Incluso si un integrador convergiera exactamente a la solución de una ecuación, dicha ecuación podría usar una deriva equivocada, una fuente de ruido irrealista o una condición de frontera inadecuada. La verificación numérica comprueba la relación entre cálculo y ecuación; la validación empírica comprueba, con criterios distintos, la relación entre ecuación y fenómeno. Confundirlas produce una falsa sensación de seguridad.

140.5. Movimiento browniano geométrico como modelo con solución exacta

El movimiento browniano geométrico es una prueba particularmente útil porque su ruido depende del estado. Consideremos la ecuación de Itô

dX_t=\mu X_t\,dt+\sigma X_t\,dW_t,\qquad X_0>0.

La solución exacta, que se obtiene aplicando la fórmula de Itô al logaritmo de $X_t$ , es

X_t=X_0\exp\!\left[\left(\mu-\frac{\sigma^2}{2}\right)t+\sigma W_t\right].

La trayectoria exacta es estrictamente positiva. Su media es $\mathbb E[X_t]=X_0e^{\mu t}$ ; la dispersión depende de $\sigma$ . El término $-\sigma^2t/2$ dentro del exponente es consecuencia directa de la variación cuadrática.

Euler–Maruyama usa la actualización

X_{n+1}^{(h)}=X_n^{(h)}[1+\mu h+\sigma\Delta W_n].

Al emplear los mismos incrementos para construir $W_T$ y la trayectoria aproximada podemos medir el error fuerte sin mezclar ruidos diferentes. En este modelo también hay una expresión exacta para la media numérica

\mathbb E[X_N^{(h)}]=X_0(1+\mu h)^N,\qquad N=T/h.

La diferencia con $X_0e^{\mu T}$ muestra un sesgo débil de orden uno, en general, para el observable $\varphi(x)=x$ . Usando $\log(1+\mu h)=\mu h-\mu^2h^2/2+O(h^3)$ , la media numérica se aproxima a la exacta con una discrepancia global proporcional a $h$ cuando $T$ permanece fijo.

La convergencia asintótica puede ocultar que la actualización de Euler–Maruyama produce $X_{n+1}<0$ , aunque $X_n>0$ , cuando $1+\mu h+\sigma\Delta W_n<0$ . Para un paso suficientemente pequeño esa posibilidad puede resultar extremadamente improbable en un horizonte determinado, pero no desaparece de la definición del método. Por eso el integrador no preserva de manera exacta la positividad que sí tiene la solución geométrica. Truncar los valores negativos a cero evitaría números no físicos a costa de introducir otra regla y, potencialmente, otro sesgo.

Una alternativa en este caso específico es usar la actualización exacta de cada intervalo

X_{n+1}=X_n\exp\!\left[\left(\mu-\frac{\sigma^2}{2}\right)h+\sigma\Delta W_n\right].

Esta expresión conserva la positividad y reproduce exactamente la solución del modelo en los instantes de la malla, porque depende de los mismos incrementos brownianos. No debemos generalizarla a difusiones no lineales sin solución analítica. Es, precisamente, el patrón contra el cual se evalúa un método aproximado. Una verificación convincente cuantifica el error y muestra cómo cambia al refinar el paso.

Cuatro curvas sobre el mismo movimiento browniano muestran una trayectoria exacta y tres discretizaciones de Euler–Maruyama. — Figura 140.2. Una trayectoria del movimiento browniano geométrico y tres aproximaciones de Euler–Maruyama construidas con incrementos anidados del mismo proceso W(t). Los detalles locales dependen de la semilla; una sola realización no establece el orden del método. Semilla 14002. Ampliar figura ↗

140.6. Estabilidad numérica y varianza falsa en Ornstein–Uhlenbeck

En sistemas con retroalimentación restauradora, una malla temporal demasiado gruesa puede alterar la distribución estacionaria o incluso hacer divergir una simulación de un proceso matemáticamente estable. Consideremos

dX_t=-\kappa X_t\,dt+\sigma\,dW_t,\qquad \kappa>0.

Es el proceso de Ornstein–Uhlenbeck ya estudiado en C131. Su media decae como $e^{-\kappa t}$ y la varianza estacionaria exacta es $V_*=\sigma^2/(2\kappa)$ . La actualización de Euler–Maruyama resulta

X_{n+1}=(1-\kappa h)X_n+\sigma\sqrt h\,Z_n.

La media numérica se multiplica por $1-\kappa h$ en cada paso. Para que las perturbaciones de la media decaigan necesitamos $|1-\kappa h|<1$ , es decir, $0<\kappa h<2$ . Si $\kappa h$ supera dos, el esquema produce una inestabilidad artificial aunque el proceso continuo sea estable para todo $\kappa>0$ . Incluso dentro del intervalo estable, la distribución estacionaria puede estar sesgada.

La varianza $V_n$ del esquema cumple la recurrencia exacta

V_{n+1}=(1-\kappa h)^2V_n+\sigma^2h.

En el límite estacionario numérico, cuando existe, despejamos

\boxed{V_h=\frac{\sigma^2}{2\kappa-\kappa^2h}.}

El denominador es menor que $2\kappa$ para $h>0$ , por lo cual Euler–Maruyama sobreestima la varianza estacionaria de esta ecuación. Con $\kappa=\sigma=1$ y $h=0{,}8$ , la varianza exacta es $0{,}5$ y la numérica teórica es $1/(2-0{,}8)=0{,}8333$ . Un histograma simulado podría parecer perfectamente gaussiano y estable, pero su anchura estaría equivocada en más de un sesenta por ciento.

El modelo tiene una transición exacta entre dos tiempos separados por $h$ ,

X_{n+1}=e^{-\kappa h}X_n+ \sigma\sqrt{\frac{1-e^{-2\kappa h}}{2\kappa}}\,Z_n.

Esta actualización conserva correctamente la distribución de transición del proceso en la malla. La comparación muestra que una malla discreta puede representar correctamente las estadísticas continuas; el sesgo procede de la aproximación de Euler–Maruyama.

Los casos con ruido aditivo y coeficientes suaves pueden tener un orden fuerte superior al $1/2$ genérico, de modo que una sola pendiente de error no describe todos los modelos. Si el objetivo científico es una varianza, una autocorrelación o una tasa de eventos, las pruebas deben medir directamente ese observable.

La varianza artificial de Euler crece al aumentar el paso, y dos campanas normales muestran densidades estacionarias de distinta anchura. — Figura 140.3. Para Ornstein–Uhlenbeck con κ=σ=1, la varianza estacionaria continua es 0,5, mientras que Euler–Maruyama produce V(h)=1/(2−h) para 0<h<2. A h=0,8 la varianza se eleva a 5/6; ambas densidades gaussianas representadas son fórmulas analíticas, no histogramas simulados. Ampliar figura ↗

140.7. Positividad, no linealidad y dominios que el método puede abandonar

Las ecuaciones usadas en ciencias de la complejidad suelen incluir variables con restricciones, como densidades de población no negativas, proporciones entre cero y uno, probabilidades normalizadas o concentraciones químicas no negativas. Una ecuación estocástica puede preservar matemáticamente esas restricciones y un esquema explícito puede violarlas. Cuando ello sucede, la primera reacción no debería ser recortar los resultados en silencio, porque ese recorte define otro algoritmo y puede cambiar las distribuciones.

Consideremos un modelo de población cuya deriva incluye una regulación fuerte y cuyo coeficiente de ruido se anula en cero. La solución exacta puede permanecer en el semieje positivo, pero una actualización explícita suficientemente grande puede saltar a valores negativos antes de que el mecanismo restaurador o la degeneración del ruido tengan el efecto esperado. En ese contexto, una trayectoria negativa no es automáticamente un hallazgo biológico; puede revelar una deficiencia de la discretización o una formulación inadecuada de las fronteras.

También hay ecuaciones con derivas superlineales para las cuales Euler–Maruyama explícito no tiene la convergencia en momentos que garantizarían las hipótesis clásicas de Lipschitz global. En esos casos, el algoritmo sin modificaciones puede producir errores numéricos extremos y resultados que empeoran al refinar el paso, aunque el proceso esté matemáticamente bien definido. Una familia de métodos denominados tamed Euler controla el crecimiento efectivo de la deriva; otros enfoques usan pasos adaptativos, esquemas implícitos o cambios de variables que preservan el dominio.

No hay una receta universal. Transformar una variable positiva mediante un logaritmo puede resolver algunos problemas, pero exige aplicar correctamente la fórmula de Itô al cambio de coordenadas. El esquema obtenido no equivale, por lo general, a simular en la variable original y tomar después el logaritmo. En sistemas con fronteras absorbentes se puede detener una trayectoria al alcanzar un umbral; en sistemas reflectantes hay que representar una dinámica compatible con la reflexión especificada, que un cambio arbitrario de signo de los valores fuera del dominio no reproduce.

El paso adaptativo es otra fuente de dificultad. Elegir el siguiente paso según el estado o el incremento aleatorio observado puede modificar inadvertidamente las propiedades de la fuente browniana si el procedimiento no está diseñado con cuidado. En problemas deterministas, el control local de error de un integrador adaptativo es familiar. En procesos estocásticos, el control debe respetar tanto la estructura probabilística como la relación entre incrementos de distintos tamaños.

Por tanto, una prueba de conservación de dominio debe formar parte de la auditoría siempre que el modelo imponga restricciones. Debemos verificar cuántas veces aparecen valores imposibles, cómo dependen de $h$ y qué ocurre cuando se usa un algoritmo que preserve matemáticamente la propiedad. La ausencia de fallos en una ejecución breve no prueba una garantía global; una simulación larga o un evento raro podrían exhibirlos después.

140.8. Primer paso, eventos raros y cruces ocultos entre puntos de la malla

En C137 y C138 usamos tiempos de primera llegada y tasas de escape. Estos observables son especialmente sensibles a cómo se observa la trayectoria entre instantes discretos. Una difusión continua puede cruzar una barrera y regresar antes de que ocurra el siguiente punto registrado. Si un programa comprueba únicamente si $X_n$ se encuentra al otro lado de la frontera, omitirá parte de esos eventos. Aumentar el número de trayectorias disminuye la incertidumbre muestral, pero no elimina el sesgo de observación entre pasos.

El problema se ve incluso en un movimiento browniano sencillo. Condicionadas a dos extremos $x$ e $y$ por debajo de un umbral $B$ , las trayectorias intermedias no están fijadas. Para un movimiento browniano con volatilidad constante $\sigma$ y sin fronteras adicionales, la probabilidad condicional de que el puente browniano atraviese $B$ en un intervalo de duración $h$ es

P\!\left(\max_{0\leq s\leq h}X_s\geq B\mid X_0=x,X_h=y\right) =\exp\!\left[-\frac{2(B-x)(B-y)}{\sigma^2h}\right],

si ambos extremos están por debajo del umbral. Esta identidad deriva del principio de reflexión para el puente browniano. En un proceso de difusión con coeficientes variables no puede usarse sin examinar sus supuestos. Muestra, sin embargo, por qué surge el problema, ya que dos observaciones que no cruzan la barrera no prueban que la trayectoria entre ellas tampoco lo haya hecho.

Una corrección basada en puentes brownianos puede reducir el sesgo para modelos adecuados. Otra opción es estudiar la convergencia de la distribución de tiempos de cruce al refinar el paso. Si la aplicación depende de cuantiles extremos, colas de supervivencia o probabilidad de llegada antes de una fecha límite, comparar solo la media es insuficiente.

Tomemos una barrera situada en $B=1$ , con una trayectoria browniana de volatilidad $\sigma=1$ que comienza en $x=0{,}2$ y termina en $y=0{,}4$ después de una unidad de tiempo. Ninguno de los extremos ha alcanzado el nivel de salida. Sin embargo, la fórmula del puente da una probabilidad condicional de cruce igual a $e^{-0{,}96}\approx0{,}383$ . Si solo se observaran los extremos, quedaría sin registrar un evento que ocurre en una proporción considerable de los puentes compatibles con ellos. Si el intervalo fuese de $h=0{,}25$ y mantuviéramos los mismos extremos, el exponente sería $-3{,}84$ , con probabilidad aproximada de $0{,}0215$ . El refinamiento disminuye el riesgo de omisión, pero el problema persiste con cualquier paso finito.

Para usar esta corrección en una simulación de primer paso, debemos establecer también qué sucede después de un cruce. Si la frontera es absorbente, la trayectoria se detiene allí y no puede contarse otra vez. Si interesa el primer cruce de una señal continua, conservar solo el signo de los estados observados puede generar errores sistemáticos en la fase de los eventos. En un problema de resonancia estocástica, esta desviación puede traducirse en un pico espurio, un desfase artificial o un aparente cambio de la intensidad óptima del ruido. Por ello conviene informar cómo se localizaron los cruces y comparar los resultados con una resolución temporal mayor, aun cuando el histograma estacionario parezca estable.

Los eventos raros añaden incertidumbre de Monte Carlo. Si un escape tiene probabilidad $p$ muy pequeña por realización, estimarlo mediante una proporción requiere un número de observaciones mucho mayor que el necesario para una media ordinaria. Para un conteo binomial de $M$ trayectorias independientes, la desviación estándar aproximada de la proporción es $\sqrt{p(1-p)/M}$ . Cuando se observan cero eventos, no podemos concluir que la probabilidad verdadera sea exactamente cero. Solo sabemos que los datos disponibles no resolvieron el fenómeno a la escala de probabilidad buscada.

Hay métodos especializados de muestreo de eventos raros, como técnicas de importancia, partición o trayectorias condicionadas. Pueden reducir costes, pero incorporan elecciones que también deben verificarse. En cualquier caso, las cifras sobre escape necesitan dos auditorías separadas, una del sesgo temporal del integrador y otra de la incertidumbre estadística de los eventos observados.

A la izquierda, varias trayectorias condicionadas con los mismos extremos bajo una barrera; algunas la cruzan. A la derecha, la probabilidad exacta de cruce disminuye con el paso. — Figura 140.4. Puentes brownianos simulados con X(0)=0,2, X(h)=0,4 y volatilidad σ=1 (izquierda, h=1), con cruces de B=1 destacados en naranja si fueron vistos en la malla interna. La curva derecha es exacta para una trayectoria continua condicionada a los extremos: P(cruce|x,y)=exp[−2(B−x)(B−y)/(σ²h)]. Resulta ≈0,383 con h=1 y ≈0,0215 con h=0,25. Semilla 14005. Ampliar figura ↗

140.9. Protocolo reproducible con las mismas trayectorias en varias resoluciones

Proponemos una verificación con dos modelos que tienen soluciones analíticas. El primero es el movimiento browniano geométrico, para el cual estimaremos el error fuerte de Euler–Maruyama con el mismo ruido a diferentes resoluciones. El segundo es Ornstein–Uhlenbeck, con una comparación entre la varianza estacionaria exacta y la varianza producida por Euler–Maruyama en un paso relativamente grande.

La secuencia comienza generando incrementos brownianos sobre una malla fina. Para obtener un incremento de una malla más gruesa, sumamos incrementos finos consecutivos; así conservamos la misma trayectoria $W_t$ en los puntos compartidos. Calculamos la solución geométrica exacta a partir de $W_T$ y repetimos el esquema de Euler–Maruyama con varios tamaños de paso. La cantidad que promediamos es $(X_T-X_N^{(h)})^2$ , cuya raíz estima el error fuerte.

import numpy as np

rng = np.random.default_rng(20261012)
T, x0, mu, sigma = 1.0, 1.0, 0.3, 0.6
paths, finest = 20000, 256
h_fine = T / finest
increments = np.sqrt(h_fine) * rng.standard_normal((paths, finest))
W_T = increments.sum(axis=1)
exact = x0 * np.exp((mu - 0.5 * sigma**2) * T + sigma * W_T)

steps = np.array([8, 16, 32, 64, 128])
strong_errors = []
for n in steps:
    h = T / int(n)
    dW = increments.reshape(paths, int(n), finest // int(n)).sum(axis=2)
    approx = x0 * np.prod(1.0 + mu*h + sigma*dW, axis=1)
    strong = np.sqrt(np.mean((approx - exact)**2))
    strong_errors.append(strong)
    weak_bias = x0 * ((1 + mu*h)**int(n) - np.exp(mu*T))
    print(f'n={n:3d}  h={h:.6f}  error fuerte={strong:.6f}'
          f'  sesgo exacto de la media={weak_bias:.6f}')

slope = np.polyfit(np.log(T / steps), np.log(strong_errors), 1)[0]
print('Pendiente fuerte estimada:', round(float(slope), 3))
assert strong_errors[-1] < strong_errors[0]

kappa, sigma_ou, h = 1.0, 1.0, 0.8
v_exact = sigma_ou**2 / (2 * kappa)
v_euler = sigma_ou**2 / (2*kappa - kappa**2*h)
print('OU: varianza exacta =', v_exact)
print('OU: varianza estacionaria de Euler =', round(v_euler, 6))
assert v_euler > v_exact

Este código verifica que el refinamiento temporal reduce el error de una trayectoria aproximada cuando se conserva el movimiento browniano de referencia. La pendiente estimada por regresión logarítmica es un indicador del orden asintótico, sin valor de demostración matemática. Para una auditoría formal de resultados numéricos convendría repetirla con más tamaños de paso, semillas independientes y errores estándar.

En la parte de Ornstein–Uhlenbeck, las dos varianzas se calculan analíticamente. No son valores muestrales; dan un criterio exacto para detectar el sesgo numérico. Si quisiéramos simular el proceso, deberíamos comprobar que las varianzas empíricas convergen a los valores previstos para cada esquema y especificar el período de calentamiento, el número de trayectorias y el manejo de las correlaciones temporales.

Para interpretar correctamente la pendiente logarítmica también debemos evitar una confusión estadística frecuente. Los puntos de las curvas de error producidas con un mismo conjunto de trayectorias acopladas son dependientes entre sí. Esa dependencia es ventajosa porque reduce la variabilidad al comparar resoluciones, pero impide tratar los puntos de la gráfica como mediciones independientes de un experimento ordinario. Si se desean bandas de incertidumbre para la pendiente, pueden repetirse lotes independientes del experimento completo y recalcular el ajuste en cada lote. Otra posibilidad es usar intervalos basados en la dispersión de los errores individuales, siempre que se conserve explícitamente el acoplamiento por trayectoria. El número de cifras significativas debe ajustarse a esa precisión, con independencia de los decimales que imprima el programa.

Una simulación publicable requiere documentar, además de la semilla, la versión del código y las bibliotecas empleadas, los parámetros, las unidades, la definición del observable y el criterio de convergencia. Es recomendable conservar un conjunto pequeño de pruebas automáticas que compare con soluciones conocidas. Cuando no hay solución exacta, la convergencia por refinamiento y la comparación con métodos independientes adquieren mayor importancia, aunque nunca sustituyen una referencia analítica cuando está disponible.

Dos gráficas logarítmicas: el error fuerte disminuye con una pendiente cercana a medio y el sesgo de la media con pendiente cercana a uno. — Figura 140.5. Error fuerte obtenido con 16 000 realizaciones acopladas de movimiento browniano geométrico y sesgo débil de la media calculado exactamente como |(1+μh)^(1/h)−e^μ| con T=1. Las pendientes ajustadas son 0,515 y 0,993. Son verificaciones para este ejemplo, no una demostración general. Ampliar figura ↗

140.10. Evaluación numérica de ecuaciones estocásticas

Los resultados de un integrador deben auditarse al nivel de cada observable. Un error pequeño en la distribución final no asegura precisión en la probabilidad de cruzar un umbral; tampoco demuestra que la trayectoria aproximada siga de cerca la realización exacta. Para estudiar ambos errores conviene reutilizar los mismos incrementos brownianos, variar sistemáticamente la malla y contrastar casos con soluciones analíticas. Una simulación puede reproducir una media y, sin embargo, tener una varianza estacionaria sesgada.

La actualización esencial del algoritmo es breve. Su interpretación no lo es. La escala $\sqrt h$ del incremento aleatorio conserva la variación cuadrática; el uso del extremo izquierdo corresponde a Itô. Cuando se comparan trayectorias, los métodos deben compartir el ruido de referencia. Cuando se comparan distribuciones, hay que separar el sesgo de discretización del error de Monte Carlo. Si se estudia un umbral, se debe reconocer que un cruce puede ocurrir entre dos instantes observados.

La distinción entre error fuerte y débil es decisiva para definir una simulación satisfactoria. Un método puede reproducir la media final de un sistema con trayectorias individuales poco precisas; las pruebas para una tasa de escape tampoco son idénticas a las pruebas para una distribución marginal. La verificación se diseña desde el objetivo científico.

Euler–Maruyama puede violar restricciones cualitativas del modelo. Positividad, conservación, estabilidad o absorción en una frontera no se heredan necesariamente de la ecuación continua. Algunos fallos se corrigen refinando el paso; otros requieren métodos especialmente construidos. El coste computacional de corregirlos forma parte de la evaluación del modelo.

La comparación de trayectorias con los mismos incrementos brownianos recuperó una pendiente de error fuerte cercana a $1/2$ en el modelo geométrico. La aproximación de la esperanza muestra un error de otro orden; el proceso de Ornstein–Uhlenbeck exhibe el sesgo de varianza que puede persistir en simulaciones aparentemente estables. Son casos concretos; en ecuaciones no lineales el comportamiento puede ser distinto.

Un informe numérico reproducible debería especificar el proceso, el observable, la convención estocástica, la malla, los incrementos utilizados y la incertidumbre de Monte Carlo. Debe comprobar también que el método respeta las restricciones cualitativas pertinentes. Cuando se investigan eventos raros, el control del sesgo de detección de cruces puede importar más que la precisión de una media global.

Referencias

Kloeden, P. E., y Platen, E. (1992). Numerical Solution of Stochastic Differential Equations. Springer. DOI: 10.1007/978-3-662-12616-5.
Higham, D. J. (2001). An algorithmic introduction to numerical simulation of stochastic differential equations. SIAM Review, 43(3), 525–546. DOI: 10.1137/S0036144500378302.
Milstein, G. N. (1995). Numerical Integration of Stochastic Differential Equations. Springer.
Øksendal, B. (2003). Stochastic Differential Equations: An Introduction with Applications (6.ª ed.). Springer.
Hutzenthaler, M., Jentzen, A., y Kloeden, P. E. (2011). Strong and weak divergence in finite time of Euler's method for stochastic differential equations with non-globally Lipschitz continuous coefficients. Proceedings of the Royal Society A, 467, 1563–1576. DOI: 10.1098/rspa.2010.0348.
Glasserman, P. (2004). Monte Carlo Methods in Financial Engineering. Springer. Simulación de movimientos brownianos y métodos de puente.