6  Métodos Numéricos para EDE

6.1 Introducción a la Discretización de EDE con difusion y salto.

Las Ecuaciones Diferenciales Estocásticas (EDE) constituyen una herramienta fundamental en la modelización matemática de fenómenos dinámicos sometidos a incertidumbre. En particular, las ecuaciones diferenciales estocásticas con saltos, por sus siglas en inglés Jump-Diffusion Stochastic Differential Equations que denotaremos de ahora en adelante por (EDSEJ) extienden el marco clásico de difusión (dominado por el movimiento browniano) al incorporar procesos de Lévy o procesos de Poisson compuestos, lo que permite capturar discontinuidades abruptas, eventos extremos o choques externos.

A pesar del sólido desarrollo teórico existente en torno a la existencia, unicidad y propiedades cualitativas de soluciones de EDSEJs, la obtención de soluciones explícitas es usual. Por ello, los métodos numéricos se vuelven indispensables tanto para la simulación como para el análisis cuantitativo de estos modelos. Sin embargo, la presencia simultánea de componentes de difusión continua y saltos discontinuos introduce desafíos computacionales y teóricos sustanciales que no están presentes en el caso puramente difusivo.

Las ecuaciones diferenciales estocásticas (SDEs) con saltos constituyen una clase fundamental de modelos en finanzas matemáticas, teoría de seguros y dinámica de poblaciones. La forma general de una SDE con saltos en \(\mathbb{R}^d\) es:

\[\begin{aligned} X_t = {}& X_0 + \int_0^t a(s,X_s)\,d s + \int_0^t b(s,X_s)\,d W_s + \int_0^t \int_E R(s,X_{s-},v)\,\tilde{p}_\phi(d v, ds) \\ & + \int_0^t \int_E R(s,X_{s-},v)\,\nu(d v)\,ds, \end{aligned} \tag{6.1}\] donde:

  • \(a: [0,T] \times \mathbb{R} \to \mathbb{R}\) es el coeficiente de deriva (drift),

  • \(b: [0,T] \times \mathbb{R} \to \mathbb{R}\) es el coeficiente de difusión,

  • \(W = \{W_t\}_{t \geq 0}\) es un proceso de Wiener, - \(p_\phi\) es una medida aleatoria de Poisson con intensidad \(\nu(dv) dt\),

  • \(\tilde{p}_\phi(dv, ds) = p_\phi(dv, ds) - \nu(dv) ds\) es la medida de Poisson compensada,

  • \(R: [0,T] \times \mathbb{R} \times E \to \mathbb{R}\) describe la magnitud de los saltos,

  • \(E \subseteq \mathbb{R}^k \setminus \{0\}\) es el espacio de marcas (marks).

La SDE con saltos compensados es: \[dX_t = \underbrace{a(t, X_t)\,dt}_{\text{drift original}} + \underbrace{b(t, X_t)\,dW_t}_{\text{difusión}} + \underbrace{\int_E R(t, X_{t-}, v)\,\tilde{p}_\phi(dv, dt)}_{\text{saltos compensados}}\] donde la medida compensada es: \[\tilde{p}_\phi(dv, dt) = p_\phi(dv, dt) - \phi(dv)\,dt\] Al expandir la integral con la medida compensada, obtenemos: \[\int_E R(t, X_{t-}, v)\,\tilde{p}_\phi(dv, dt) = \underbrace{\int_E R(t, X_{t-}, v)\,p_\phi(dv, dt)}_{\text{saltos brutos}} - \underbrace{\int_E R(t, X_{t-}, v)\,\phi(dv)\,dt}_{\text{compensación}}\] Por lo tanto, la EDE completa se puede escribir como: \[dX_t = a(t, X_t)\,dt + b(t, X_t)\,dW_t + \int_E R(t, X_{t-}, v)\,p_\phi(dv, dt) - \int_E R(t, X_{t-}, v)\,\phi(dv)\,dt\] Para discretizar esta ecuación, reorganizamos los términos de la siguiente manera: \[dX_t = \underbrace{\left[a(t, X_t) - \int_E R(t, X_{t-}, v)\,\phi(dv)\right]}_{\text{drift compensado}} dt + b(t, X_t)\,dW_t + \int_E R(t, X_{t-}, v)\,p_\phi(dv, dt)\] El término de compensación cambia de signo al moverse al drift, porque: \[a(t, X_t)\,dt - \int_E R(t, X_{t-}, v)\,\phi(dv)\,dt = \left[a(t, X_t) - \int_E c(t, X_{t-}, v)\,\phi(dv)\right] dt\]

6.1.1 Formulación general del modelo de Cramér–Lundberg con reinversión como EDESJ.

El modelo de Cramér–Lundberg extendido con reinversión y parámetros temporales puede escribirse como una ecuación diferencial estocástica con saltos (EDE). Consideremos la siguiente ecuación diferencial estocástica con saltos (EDE): \[d X_t = \underbrace{c(t)\,dt}_{\text{primas}} + \underbrace{r(t) X_t\,dt}_{\text{rendimiento}} + \underbrace{\sigma(t) X_t\,d W_t}_{\text{volatilidad}} - \underbrace{\int_E R(t, X_{t-}, v)\,\tilde{p}_\phi(dv, dt)}_{\text{saltos compensados}}, \tag{6.2}\] donde

  • \(c(t) \in \mathbb{R}\) es un flujo de primas determinista y acotado,

  • \(r(t), \sigma(t)\) son funciones acotadas y medibles,

  • \(R(t, x, v)\) es la magnitud del salto, dependiente del tiempo \(t\), del estado \(x\) y de la marca \(v\), - \(\tilde{p}_\phi(dv, d t) = p_\phi(d v, dt) - \phi(dv)d t\) es la medida de Poisson compensada, - \(\phi(E) < \infty\) y \(\int_E |R(t, x, v)|^2 \phi(d v) < \infty\) para todo \(t \in [0,T]\), \(x \in \mathbb{R}\).

6.2 Esquema de Milstein

Sea la ecuación diferencial estocástica con saltos: \[dX_t = a(t, X_t)\,dt + b(t, X_t)\,dW_t + \int_E R(t, X_{t-}, v)\, p_\phi(dv, dt), \tag{1}\] donde \(p_\phi\) es la medida aleatoria de Poisson con intensidad \(\phi(dv)dt\).

La discretización de Milstein sin compensación se define como: \[\begin{aligned} Y_{n+1} = Y_n &+ a(t_n, Y_n),\Delta t \\ &+ b(t_n, Y_n)\,\Delta W_n \\ &+ \frac{1}{2}b(t_n, Y_n)b_x(t_n, Y_n)\left[(\Delta W_n)^2 - \Delta t\right] \\ &+ \sum_{i=1}^{\Delta N_n} R(t_n, Y_n, \xi_i), \end{aligned}\] con:

  • \(\Delta t = t_{n+1} - t_n\) paso de tiempo,

  • \(\Delta W_n \sim \mathcal{N}(0, \Delta t)\) incremento browniano,

  • \(\Delta N_n \sim \text{Poisson}(\lambda \Delta t)\), \(\lambda = \phi(E)\), número de saltos en el intervalo,

  • \(\xi_i \overset{\text{i.i.d.}}{\sim} \phi(dv)/\lambda\) marcas de los saltos,

  • \(b_x = \partial b/\partial x\) derivada del coeficiente de difusión,

  • \(a(t,x)\) y \(b(t,x)\) coeficientes de deriva y difusión.

Teorema 6.1 (Convergencia fuerte de orden \(\gamma = 1/2\)). Bajo las hipótesis:

  1. Lipschitz global: \(\|a(t, x_1) - a(t, x_2)\| \leq L_a \|x_1 - x_2\|\), análogo para \(b\) y \(R\).
  2. Crecimiento lineal: \(\|a\| + \|b\| + \left(\int_E \|R\|^2 \phi\right)^{1/2} \leq K(1 + \|x\|)\).
  3. Momentos finitos de saltos: \(\int_E \|R(t, x, v)\|^4 \phi(dv) < \infty\).
  4. Condiciones iniciales: \(\mathbb{E}[\|X_0\|^2] < \infty\) y \(Y_0 = X_0\).

Entonces, existe una constante \(C > 0\), independiente de \(\Delta t\), tal que:

\[ \max_{0 \leq n \leq N} \left(\mathbb{E}[\|X_{t_n} - Y_n\|^2]\right)^{1/2} \leq C \sqrt{\Delta t}. \tag{6.3}\]

Es decir, el esquema converge fuertemente con orden \(\gamma = 1/2\).

6.3 Esquema de Milstein Compensado

Sea \(Y_n\) la aproximación numérica dada por el esquema Milstein compensado: \[\begin{aligned} Y_{n+1} = Y_n &+ \left[ a(t_n, Y_n) - \int_E R(t_n, Y_n, v) \phi(d v) \right]\Delta t \\ &+ b(t_n, Y_n) \Delta W_n \\ &+ \frac{1}{2} b(t_n, Y_n) b_x(t_n, Y_n) \left[ (\Delta W_n)^2 - \Delta t \right] \\ &+ \sum_{i=1}^{\Delta N_n} R(t_n, Y_n, \xi_i), \end{aligned} \tag{6.4}\]

  • \(\Delta t = t_{n+1} - t_n\) es el paso de tiempo,
  • \(\Delta W_n \sim \mathcal{N}(0, \Delta t)\) es el incremento del proceso de Wiener,
  • \(\Delta N_n \sim \text{Poisson}(\lambda \Delta t)\), \(\lambda = \phi(E)\), es el número de saltos en el intervalo,
  • \(\xi_i \stackrel{\text{i.i.d.}}{\sim} \phi(dv)/\lambda\) son las marcas de los saltos,
  • \(b_x = \partial b / \partial x\) es la derivada parcial del coeficiente de difusión,
  • \(a(t, x) = c(t) + r(t)x\) y \(b(t, x) = \sigma(t)x\).

(Convergencia fuerte de orden \(\gamma = 1/2\) ). Supongamos las siguientes hipótesis:

  1. Lipschitz global: Los coeficientes \(a(t,x)\) y \(b(t,x)\) son globalmente Lipschitz en \(x\), y \(c(t,x,v)\) es globalmente Lipschitz en \(x\) uniformemente en \(t\) y \(v\): \[|c(t, x_1, v) - c(t, x_2, v)| \leq L_c |x_1 - x_2|, \quad \forall t \in [0,T], \, x_1, x_2 \in \mathbb{R}, \, v \in E.\]

  2. Crecimiento lineal: Existe \(K > 0\) tal que \[|a(t,x)| + |b(t,x)| + \left( \int_E |R(t, x, v)|^2 \phi(d v) \right)^{1/2} \leq K(1 + |x|),\] para todo \(t \in [0,T]\) y \(x \in \mathbb{R}\).

  3. Momentos finitos de saltos: \(\int_E |R(t, x, v)|^4 \phi(d v) < \infty\) para todo \(t \in [0,T]\) y \(x \in \mathbb{R}\).

  4. Condiciones iniciales: \(\mathbb{E}[|X_0|^2] < \infty\) y \(Y_0 = X_0\). Entonces, existe una constante \(C > 0\), independiente de \(\, dt\), tal que: \[\max_{0 \leq n \leq N} \left( \mathbb{E}\left[ |X_{t_n} - Y_n|^2 \right] \right)^{1/2} \leq C \sqrt{dt}. \tag{6.5}\] Es decir, el esquema converge fuertemente con orden \(\gamma = 1/2\).

Demostración. La EDE Ecuación 6.2 puede reescribirse en la forma estándar Platen y Bruti-Liberati (2010), Ec. (1.8.2): \[d X_t = a(t,X_t)\,d t + b(t,X_t)\,dW_t + \int_E R(t,X_{t-},v)\,p_\phi(dv, d t), \tag{6.6}\]

donde:

  • \(a(t,x) = c(t) + r(t)x\) (drift original),

  • \(b(t,x) = \sigma(t)x\),

  • \(R(t,x,v)\) es el tamaño del salto bruto, dependiente del estado. Las hipótesis de Lipschitz y crecimiento lineal garantizan que los coeficientes satisfacen las condiciones del Teorema 1.9.3 de Platen y Bruti-Liberati (2010), que establece la existencia y unicidad de una solución fuerte \(X_t\).

Además, el Teorema 1.9.3 garantiza que la solución tiene momentos finitos de segundo orden: \[\mathbb{E}\left[ \sup_{0 \leq t \leq T} |X_t|^2 \right] \leq C_1 (1 + \mathbb{E}[|X_0|^2]) < \infty. \tag{6.7}\]

El esquema Ecuación 6.4 incluye la corrección de Milstein, que corresponde a truncar la expansión estocástica de Wagner-Platen (ver Platen y Bruti-Liberati (2010), Ec. (4.4.7)) después de los términos de orden \(\mathcal{O}((dt)^{3/2})\).

Para una SDE con saltos, la expansión completa incluye términos adicionales relacionados con la medida de Poisson. Como el modelo es aditivo en los saltos (aunque \(c(t,x,v)\) depende de \(x\), no hay términos cruzados de orden superior), los operadores de salto y difusión conmutan trivialmente, lo que simplifica el análisis.

Definimos el error local \(R_n\) como la diferencia entre la solución exacta de la EDE avanzada un paso de tiempo \(\Delta t\) y el valor producido por el esquema numérico, ambos partiendo del mismo punto inicial \(Y_n\) en el tiempo \(t_n\):

\[ R_n := X_{t_{n+1}}^{t_n, Y_n} - Y_{n+1}, \tag{6.8}\]

donde \(X_{t_{n+1}}^{t_n, Y_n}\) denota la solución de la EDE Ecuación 6.2 evaluada en \(t_{n+1}\), con condición inicial \(X_{t_n} = Y_n\).

El esquema Ecuación 6.4 se reescribe explícitamente como:

\[ Y_{n+1} = Y_n + a(t_n, Y_n)\Delta t + b(t_n, Y_n)\Delta W_n + \frac{1}{2}b(t_n, Y_n)b'(t_n, Y_n)[(\Delta W_n)^2 - \Delta t] + \int_E R(t_n, Y_n, v)\tilde{p}_\phi(dv, \Delta t) + C_n \Delta t, \tag{6.9}\]

donde:

  • \(\Delta W_n = W_{t_{n+1}} - W_{t_n} \sim \mathcal{N}(0, \Delta t)\) es el incremento del movimiento browniano,

  • \(\tilde{p}_\phi(dv, \Delta t) = p_\phi(dv, [t_n, t_{n+1}]) - \phi(dv)\Delta t\) es la medida de Poisson compensada,

  • \(C_n = -\int_E R(t_n, Y_n, v)\phi(dv) = -\lambda \mathbb{E}[Y]\) es el término de compensación que garantiza \(\mathbb{E}\left[\int_E R \tilde{p}_\phi\right] = 0\).

Por el teorema de expansión estocástica de Wagner-Platen (platen2010?, Ecuación (4.4.7)), la solución exacta \(X_{t_{n+1}}^{t_n, Y_n}\) admite la expansión:

\[ X_{t_{n+1}}^{t_n, Y_n} = Y_n + a\Delta t + b\Delta W_n + \frac{1}{2}bb'[(\Delta W_n)^2 - \Delta t] + \int_E R\tilde{p}_\phi + \rho_n, \tag{6.10}\]

donde todos los coeficientes \(a, b, b', R\) se evalúan en \((t_n, Y_n)\), y \(\rho_n\) es el residuo de orden superior que contiene todos los términos de la expansión de Itô-Taylor de orden \((\Delta t)^{3/2}\) o superior. Explícitamente:

\[\begin{aligned} \rho_n = \underbrace{\int_{t_n}^{t_{n+1}} \int_{t_n}^{s} \mathcal{L}^0 a \, du \, ds}_{\text{término de deriva de orden } (\Delta t)^2} &+ \underbrace{\int_{t_n}^{t_{n+1}} \int_{t_n}^{s} \mathcal{L}^1 b \, dW_u \, ds}_{\text{término cruzado difusión-deriva}} \\ &+ \underbrace{\int_{t_n}^{t_{n+1}} \int_E [\mathcal{L}^{-1} R] \tilde{p}_\phi(dv, ds)}_{\text{término de salto de orden superior}} \\ &+ \cdots \end{aligned} \tag{6.11}\]

Donde \(\mathcal{L}^0\) y \(\mathcal{L}^1\) son los operadores diferenciales de Itô definidos por:

\[ \mathcal{L}^0 := \frac{\partial}{\partial t} + a(t,x)\frac{\partial}{\partial x} + \frac{1}{2}b(t,x)^2\frac{\partial^2}{\partial x^2}, \tag{6.12}\]

\[ \mathcal{L}^1 := b(t,x)\frac{\partial}{\partial x}. \tag{6.13}\]

Estos operadores surgen de aplicar la fórmula de Itô iteradamente a los coeficientes de la EDE. El operador \(\mathcal{L}^0\) actúa sobre las funciones deterministas y produce términos de orden \((\Delta t)^2\), mientras que \(\mathcal{L}^1\) actúa sobre las funciones estocásticas y produce términos de orden \((\Delta t)^{3/2}\) cuando se integra contra \(dW\).

Al restar Ecuación 6.9 de Ecuación 6.10, todos los términos de orden \(\Delta t\) y \((\Delta W_n)^2 - \Delta t\) se cancelan exactamente, quedando:

\[ R_n = X_{t_{n+1}}^{t_n, Y_n} - Y_{n+1} = \rho_n. \tag{6.14}\]

Por tanto, el error local coincide exactamente con el residuo de la expansión de Itô-Taylor.

Para demostrar que \(\mathbb{E}[\|R_n\|^2 \mid \mathcal{F}_{t_n}] \leq K(1 + \|Y_n\|^2)(\Delta t)^2\). Por Ecuación 6.14, basta acotar \(\mathbb{E}[\|\rho_n\|^2 \mid \mathcal{F}_{t_n}]\).

Analizando cada termino Ecuación 6.11 por separaddo y aplicando las desigualdades de isometría estocástica se tiene:

\[ I_1 := \int_{t_n}^{t_{n+1}} \int_{t_n}^{s} \mathcal{L}^0 a \, du \, ds. \] Bajo las hipótesis de Lipschitz y crecimiento lineal, \(\mathcal{L}^0 a\) está acotado por \(K(1 + \|Y_n\|)\). Aplicando la desigualdad de Cauchy-Schwarz para integrales deterministas:

\[\begin{aligned} \mathbb{E}[|I_1|^2 \mid \mathcal{F}_{t_n}] \leq \left(\int_{t_n}^{t_{n+1}} \int_{t_n}^{s} K(1 + \|Y_n\|) \, du \, ds\right)^2 = K^2(1 + \|Y_n\|^2)\left(\frac{(\Delta t)^2}{2}\right)^2 &= K_1(1 + \|Y_n\|^2)(\Delta t)^4. \end{aligned} \tag{6.15}\]

\[ I_2 := \int_{t_n}^{t_{n+1}} \int_{t_n}^{s} \mathcal{L}^1 b \, dW_u \, ds. \] Aplicaando isometría de Itô (protter2005?, Teorema IV.28): para una integral estocástica \(\int_0^t H_u \, dW_u\) con \(H\) adaptado y \(\mathbb{E}\int_0^t |H_u|^2 du < \infty\), se tiene \(\mathbb{E}\left[\left|\int_0^t H_u dW_u\right|^2\right] = \mathbb{E}\int_0^t |H_u|^2 du\). Aplicando esto condicionalmente a \(\mathcal{F}_{t_n}\): \[ \mathbb{E}[|I_2|^2 \mid \mathcal{F}_{t_n}] = \mathbb{E}\left[\int_{t_n}^{t_{n+1}} \left|\int_{t_n}^{s} \mathcal{L}^1 b \, dW_u\right|^2 ds \;\middle|\; \mathcal{F}_{t_n}\right]. \] Usando nuevamente la isometría de Itô en la integral interior: \[ \mathbb{E}[|I_2|^2 \mid \mathcal{F}_{t_n}] = \int_{t_n}^{t_{n+1}} \int_{t_n}^{s} \mathbb{E}[|\mathcal{L}^1 b|^2 \mid \mathcal{F}_{t_n}] \, du \, ds \leq K^2(1 + \|Y_n\|^2) \int_{t_n}^{t_{n+1}} (s - t_n) \, ds = K_2(1 + \|Y_n\|^2)(\Delta t)^2. \tag{6.16}\]

\[ I_3 := \int_{t_n}^{t_{n+1}} \int_E [\mathcal{L}^{-1} R] \tilde{p}_\phi(dv, ds). \] Luego aplicamos la isometría para medidas de Poisson compensadas (platen2010?, Ecuación (4.2.1)): para una integral \(\int_0^t \int_E H(v,s) \tilde{p}_\phi(dv, ds)\) con \(H\) predecible y \(\mathbb{E}\int_0^t \int_E |H|^2 \phi(dv) ds < \infty\), se tiene: \[ \mathbb{E}\left[\left|\int_0^t \int_E H \tilde{p}_\phi\right|^2\right] = \mathbb{E}\int_0^t \int_E |H|^2 \phi(dv) ds. \] Aplicando esto condicionalmente: \[ \mathbb{E}[|I_3|^2 \mid \mathcal{F}_{t_n}] = \int_{t_n}^{t_{n+1}} \int_E \mathbb{E}[|\mathcal{L}^{-1} R|^2 \mid \mathcal{F}_{t_n}] \phi(dv) ds \leq K^2(1 + \|Y_n\|^2) \Delta t \int_E \phi(dv) = K_3(1 + \|Y_n\|^2)(\Delta t)^2, \tag{6.17}\] donde usamos la hipótesis de momentos finitos de saltos (Hipótesis 3) para garantizar que \(\int_E \phi(dv) = \lambda < \infty\).

Usando la desigualdad \((a + b + c)^2 \leq 3(a^2 + b^2 + c^2)\) y sumando Ecuación 6.15, Ecuación 6.16, Ecuación 6.17:

\[ \mathbb{E}[\|\rho_n\|^2 \mid \mathcal{F}_{t_n}] \leq 3\left[K_1(1 + \|Y_n\|^2)(\Delta t)^4 + K_2(1 + \|Y_n\|^2)(\Delta t)^2 + K_3(1 + \|Y_n\|^2)(\Delta t)^2\right]. \]

Para \(\Delta t \leq 1\), \((\Delta t)^4 \leq (\Delta t)^2\), luego:

\[ \mathbb{E}[\|R_n\|^2 \mid \mathcal{F}_{t_n}] = \mathbb{E}[\|\rho_n\|^2 \mid \mathcal{F}_{t_n}] \leq K(1 + \|Y_n\|^2)(\Delta t)^2, \tag{6.18}\]

donde \(K = 3(K_1 + K_2 + K_3)\). Lo cual verifica la condición de Cota del error local en \(L^2\) del Teorema Fundamental con \(\gamma = 1/2\).

Para propagar el error local al error global, necesitamos descomponer \(R_n\) en una parte determinista (sesgo) y una parte martingala (ruido). Esta descomposición es fundamental porque la parte martingala tiene esperanza condicional cero y no contribuye al sesgo acumulado.

Definimos la proyección condicional del error local:

\[ \bar{R}_n := \mathbb{E}[R_n \mid \mathcal{F}_{t_n}], \tag{6.19}\]

y la parte martingala:

\[ M_n := R_n - \bar{R}_n = R_n - \mathbb{E}[R_n \mid \mathcal{F}_{t_n}]. \tag{6.20}\]

Por construcción, \(M_n\) satisface: 1. \(\mathbb{E}[M_n \mid \mathcal{F}_{t_n}] = 0\). 2. \(R_n = \bar{R}_n + M_n\).

Usando la desigualdad de Jensen condicional \(\|\mathbb{E}[X \mid \mathcal{F}]\|^2 \leq \mathbb{E}[\|X\|^2 \mid \mathcal{F}]\):

\[ \mathbb{E}[\|M_n\|^2 \mid \mathcal{F}_{t_n}] = \mathbb{E}[\|R_n - \bar{R}_n\|^2 \mid \mathcal{F}_{t_n}] = \mathbb{E}[\|R_n\|^2 \mid \mathcal{F}_{t_n}] - \|\bar{R}_n\|^2 \leq \mathbb{E}[\|R_n\|^2 \mid \mathcal{F}_{t_n}] \leq K(1 + \|Y_n\|^2)(\Delta t)^2. \tag{6.21}\]

El sesgo \(\bar{R}_n\) proviene principalmente del término \(I_1\) en Ecuación 6.11, que es de orden \((\Delta t)^2\) en valor esperado. Los términos \(I_2\) e \(I_3\) tienen valor esperado condicional cero por ser integrales estocásticas contra \(dW\) y \(\tilde{p}_\phi\) respectivamente. Por tanto:

\[ \|\bar{R}_n\| = \|\mathbb{E}[R_n \mid \mathcal{F}_{t_n}]\| \leq K(1 + \|Y_n\|)(\Delta t)^2. \tag{6.22}\]

Definimos el error global en el paso \(n\) como:

\[ e_n := X_{t_n} - Y_n, \tag{6.23}\]

donde \(X_{t_n}\) es la solución exacta de la EDE (Ecuación 6.2) con condición inicial \(X_0\), y \(Y_n\) es la aproximación numérica dada por el esquema Ecuación 6.9.

Restando Ecuación 6.9 de la EDE integrada exactamente desde \(t_n\) hasta \(t_{n+1}\):

\[ X_{t_{n+1}} = X_{t_n} + \int_{t_n}^{t_{n+1}} a(s, X_s) ds + \int_{t_n}^{t_{n+1}} b(s, X_s) dW_s + \int_{t_n}^{t_{n+1}} \int_E R(s, X_{s^-}, v) \tilde{p}_\phi(dv, ds), \]

y usando la definición de error local Ecuación 6.8 con \(Y_n\) como punto inicial, obtenemos:

\[ e_{n+1} = e_n + \int_{t_n}^{t_{n+1}} [a(s, X_s) - a(t_n, Y_n)] ds + \int_{t_n}^{t_{n+1}} [b(s, X_s) - b(t_n, Y_n)] dW_s + R_n. \tag{6.24}\]

Para simplificar, aproximamos \(a(s, X_s) \approx a(t_n, X_{t_n})\) y \(b(s, X_s) \approx b(t_n, X_{t_n})\) para \(s \in [t_n, t_{n+1}]\) y se tiene que

\[ e_{n+1} \approx e_n + [a(t_n, X_{t_n}) - a(t_n, Y_n)]\Delta t + [b(t_n, X_{t_n}) - b(t_n, Y_n)]\Delta W_n + R_n. \tag{6.25}\]

Tomamos la norma al cuadrado de Ecuación 6.25 y aplicamos la esperanza condicional \(\mathbb{E}[\cdot \mid \mathcal{F}_{t_n}]\):

\[ \mathbb{E}[\|e_{n+1}\|^2 \mid \mathcal{F}_{t_n}] = \mathbb{E}\left[\|e_n + \Delta a_n \Delta t + \Delta b_n \Delta W_n + R_n\|^2 \;\middle|\; \mathcal{F}_{t_n}\right], \]

donde \(\Delta a_n := a(t_n, X_{t_n}) - a(t_n, Y_n)\) y \(\Delta b_n := b(t_n, X_{t_n}) - b(t_n, Y_n)\).

Expandiendo el cuadrado y usando las propiedades de independencia condicional: - \(\mathbb{E}[\Delta W_n \mid \mathcal{F}_{t_n}] = 0\) - \(\mathbb{E}[R_n \mid \mathcal{F}_{t_n}] = \bar{R}_n\) - \(\mathbb{E}[\Delta W_n \cdot R_n \mid \mathcal{F}_{t_n}] = \mathbb{E}[\Delta W_n \cdot M_n \mid \mathcal{F}_{t_n}]\),

obtenemos:

\[ \mathbb{E}[\|e_{n+1}\|^2 \mid \mathcal{F}_{t_n}] = \|e_n\|^2 + 2e_n \cdot \Delta a_n \Delta t + 2e_n \cdot \bar{R}_n + \|\Delta b_n\|^2 \Delta t + \mathbb{E}[\|R_n\|^2 \mid \mathcal{F}_{t_n}] + 2\mathbb{E}[\Delta b_n \Delta W_n \cdot M_n \mid \mathcal{F}_{t_n}]. \tag{6.26}\]

El último término se anula porque \(\Delta W_n\) y \(M_n\) son independientes condicionalmente a \(\mathcal{F}_{t_n}\).

Usando las desigualdades de:

  1. Desigualdad de Cauchy-Schwarz: \(2e_n \cdot \Delta a_n \Delta t \leq 2\|e_n\| \|\Delta a_n\| \Delta t\),
  2. Desigualdad de Young: \(2ab \leq \epsilon a^2 + \frac{1}{\epsilon} b^2\) para todo \(\epsilon > 0\),
  3. Condición de Lipschitz (Hipótesis 1): \(\|\Delta a_n\| \leq L_a \|e_n\|\) y \(\|\Delta b_n\| \leq L_b \|e_n\|\).

Aplicando estas desigualdades con \(\epsilon = 1\):

\[ 2e_n \cdot \Delta a_n \Delta t \leq \|e_n\|^2 \Delta t + L_a^2 \|e_n\|^2 \Delta t = (1 + L_a^2)\|e_n\|^2 \Delta t, \]

\[ \|\Delta b_n\|^2 \Delta t \leq L_b^2 \|e_n\|^2 \Delta t, \]

\[ 2e_n \cdot \bar{R}_n \leq \|e_n\|^2 + \|\bar{R}_n\|^2 \leq \|e_n\|^2 + K^2(1 + \|Y_n\|^2)(\Delta t)^4. \]

Sustituyendo en Ecuación 6.26 y usando la cota Ecuación 6.18:

\[ \mathbb{E}[\|e_{n+1}\|^2 \mid \mathcal{F}_{t_n}] \leq \|e_n\|^2 + (1 + L_a^2 + L_b^2)\|e_n\|^2 \Delta t + \|e_n\|^2 + K(1 + \|Y_n\|^2)(\Delta t)^2 + K^2(1 + \|Y_n\|^2)(\Delta t)^4. \]

Para \(\Delta t \leq 1\), \((\Delta t)^4 \leq (\Delta t)^2\), luego agrupando términos:

\[ \mathbb{E}[\|e_{n+1}\|^2 \mid \mathcal{F}_{t_n}] \leq (1 + K_1 \Delta t)\|e_n\|^2 + K_2(1 + \|Y_n\|^2)(\Delta t)^2, \tag{6.27}\]

donde \(K_1 = 2 + L_a^2 + L_b^2\) y \(K_2 = K + K^2\).

Tomamos el valor esperado total en Ecuación 6.27, usando la propiedad de la torre \(\mathbb{E}[\mathbb{E}[X \mid \mathcal{F}]] = \mathbb{E}[X]\):

\[ \mathbb{E}[\|e_{n+1}\|^2] \leq (1 + K_1 \Delta t)\mathbb{E}[\|e_n\|^2] + K_2(1 + \mathbb{E}[\|Y_n\|^2])(\Delta t)^2. \tag{6.28}\]

Definimos \(E_n := \mathbb{E}[\|e_n\|^2]\) y usamos la cota de momentos finitos uniformes \(\mathbb{E}[\|Y_n\|^2] \leq C_Y\) para todo \(n\). Entonces:

\[ E_{n+1} \leq (1 + K_1 \Delta t)E_n + K_3 (\Delta t)^2, \tag{6.29}\]

donde \(K_3 = K_2(1 + C_Y)\).

Luego aplicando Ecuación 6.29 recursivamente desde \(n = 0\) hasta \(n = N-1\), usando la condición inicial \(E_0 = \mathbb{E}[\|e_0\|^2] = \mathbb{E}[\|X_0 - Y_0\|^2] = 0\) tenemos:

\[ E_1 \leq (1 + K_1 \Delta t) \cdot 0 + K_3 (\Delta t)^2 = K_3 (\Delta t)^2, \]

\[ E_2 \leq (1 + K_1 \Delta t)K_3 (\Delta t)^2 + K_3 (\Delta t)^2 = K_3 (\Delta t)^2 [1 + (1 + K_1 \Delta t)], \]

\[ E_3 \leq (1 + K_1 \Delta t)K_3 (\Delta t)^2 [1 + (1 + K_1 \Delta t)] + K_3 (\Delta t)^2 = K_3 (\Delta t)^2 [1 + (1 + K_1 \Delta t) + (1 + K_1 \Delta t)^2]. \]

Por inducción, para todo \(n\):

\[ E_n \leq K_3 (\Delta t)^2 \sum_{k=0}^{n-1} (1 + K_1 \Delta t)^k. \tag{6.30}\]

y la suma geométrica se evalúa explícitamente:

\[ \sum_{k=0}^{n-1} (1 + K_1 \Delta t)^k = \frac{(1 + K_1 \Delta t)^n - 1}{(1 + K_1 \Delta t) - 1} = \frac{(1 + K_1 \Delta t)^n - 1}{K_1 \Delta t}. \]

Sustituyendo en Ecuación 6.30 con \(n = N\) y \(N \Delta t = T\):

\[ E_N \leq K_3 (\Delta t)^2 \cdot \frac{(1 + K_1 \Delta t)^N - 1}{K_1 \Delta t} = \frac{K_3}{K_1} \Delta t \left[(1 + K_1 \Delta t)^N - 1\right]. \tag{6.31}\]

Usamos la desigualdad fundamental \((1 + x)^n \leq e^{nx}\) para \(x \geq 0\):

\[ (1 + K_1 \Delta t)^N \leq e^{K_1 N \Delta t} = e^{K_1 T}. \]

Por tanto:

\[ E_N \leq \frac{K_3}{K_1} (e^{K_1 T} - 1) \Delta t = C \Delta t, \tag{6.32}\]

donde \(C = \frac{K_3}{K_1} (e^{K_1 T} - 1)\) es una constante independiente de \(\Delta t\).

Tomando la raíz cuadrada en ambos lados:

\[ \left(\mathbb{E}[\|e_N\|^2]\right)^{1/2} \leq \sqrt{C} \sqrt{\Delta t}, \]

lo que demuestra la convergencia fuerte de orden \(\gamma = 1/2\).

Ahora para verificar la condición de momentos finitos uniformes, definimos:

\[ M_n := \mathbb{E}\left[\sup_{0 \leq k \leq n} \|Y_k\|^2\right]. \tag{6.33}\]

Tomando la norma al cuadrado en Ecuación 6.9 y aplicando la desigualdad \((a + b + c + d + e + f)^2 \leq 6(a^2 + b^2 + c^2 + d^2 + e^2 + f^2)\):

\[ \|Y_{n+1}\|^2 \leq 6\|Y_n\|^2 + 6\|a(t_n, Y_n)\|^2 (\Delta t)^2 + 6\|b(t_n, Y_n)\|^2 \|\Delta W_n\|^2 + 6\left\|\frac{1}{2}bb'[(\Delta W_n)^2 - \Delta t]\right\|^2 + 6\left\|\int_E R\tilde{p}_\phi\right\|^2 + 6\|C_n\|^2 (\Delta t)^2. \tag{6.34}\]

Ahota tomamos el supremo sobre \(k \in \{0, \dots, n+1\}\) y luego la esperanza. Para los términos estocásticos, aplicamos las desigualdades de Burkholder-Davis-Gundy (protter2005?, Teorema IV.42):

  • Para la integral browniana: \(\mathbb{E}\left[\sup_{0 \leq k \leq n} \left\|\sum_{j=0}^{k} b_j \Delta W_j\right\|^2\right] \leq C_{BDG} \mathbb{E}\left[\sum_{j=0}^{n} \|b_j\|^2 \Delta t\right]\),

  • Para la integral de Poisson compensada: \(\mathbb{E}\left[\sup_{0 \leq k \leq n} \left\|\sum_{j=0}^{k} \int_E R_j \tilde{p}_\phi\right\|^2\right] \leq C_{BDG} \mathbb{E}\left[\sum_{j=0}^{n} \int_E \|R_j\|^2 \phi(dv) \Delta t\right]\).

Usando la condición de crecimiento lineal, \(\|a\|^2 + \|b\|^2 + \int_E \|R\|^2 \phi \leq K^2(1 + \|Y_n\|^2)\), se obtiene:

\[ M_{n+1} \leq 6M_n + 6K^2(1 + M_n)(\Delta t)^2 + 6C_{BDG} K^2(1 + M_n)\Delta t + 6K^2(1 + M_n)\Delta t. \]

Agrupando términos y usando \(\Delta t \leq 1\):

\[ M_{n+1} \leq (1 + K_5 \Delta t)M_n + K_6 \Delta t, \tag{6.35}\]

donde \(K_5 = 6 + 6K^2 + 6C_{BDG}K^2 + 6K^2\) y \(K_6 = 6K^2 + 6C_{BDG}K^2 + 6K^2\).

Aplicando Ecuación 6.35 recursivamente desde \(n = 0\):

\[ M_N \leq M_0 (1 + K_5 \Delta t)^N + K_6 \Delta t \sum_{k=0}^{N-1} (1 + K_5 \Delta t)^k. \]

Finalmente usando \((1 + K_5 \Delta t)^N \leq e^{K_5 T}\) y la fórmula de la suma geométrica:

\[ M_N \leq M_0 e^{K_5 T} + \frac{K_6}{K_5}(e^{K_5 T} - 1) < \infty. \tag{6.36}\]

Por tanto, el esquema tiene momentos finitos uniformes. \(\quad \blacksquare\)

6.4 Esquema Semi-Implícito Compensado

El método semi-implícito compensado se define aplicando una regla de Euler implícita en el drift y explícita en la difusión y el salto compensado \[U_{n+1} = U_n + [c_{n+1} + r_{n+1}U_{n+1}]\Delta t + \sigma_n U_n \Delta W_n - \tilde{J}_n, \tag{6.37}\]

donde \(c_{n+1} = c(t_{n+1})\), \(r_{n+1} = r(t_{n+1})\), \(\sigma_n = \sigma(t_n)\), y \(\Delta W_n = W_{t_{n+1}} - W_{t_n}\).

Reordenando algebraicamente para despejar \(U_{n+1}\):

\[U_{n+1}(1 - r_{n+1}\Delta t) = U_n(1 + \sigma_n \Delta W_n) + c_{n+1}\Delta t + \mathbb{E}[\Delta J_n] - \Delta J_n,\] lo que da lugar a la fórmula cerrada central:

\[U_{n+1} = \frac{U_n (1 + \sigma(t_n) \Delta W_n ) + c(t_n)\Delta t + \int_E c(t_n, U_n, v) \phi(dv) \Delta t - \sum_{i=1}^{\Delta N_n} c(t_n, U_n, \xi_i) }{1 - r(t_{n+1}) \Delta t} \tag{6.38}\] donde:

  • \(\Delta t = t_{n+1} - t_n\),
  • \(\Delta W_n \sim \mathcal{N}(0, \Delta t)\),
  • \(\Delta N_n \sim \text{Poisson}(\lambda \Delta t)\), \(\lambda = \phi(E)\),
  • \(\xi_i \stackrel{\text{i.i.d.}}{\sim} \phi(dv)/\lambda\).

Teorema 6.2 (Convergencia fuerte de orden \(\gamma = 1/2\)) Supongamos las siguientes hipótesis:

  • (H1) Lipschitz global: Los coeficientes \(a(t,u) = c(t) + r(t)u\) y \(b(t,u) = \sigma(t)u\) son globalmente Lipschitz en \(u\), y \(c(t,u,v)\) es globalmente Lipschitz en \(u\) uniformemente en \(t\) y \(v\): \[ |c(t, u_1, v) - c(t, u_2, v)| \leq L_c |u_1 - u_2|, \quad \forall t \in [0,T], \, u_1, u_2 \in \mathbb{R}, \, v \in E. \tag{H1} \]

  • (H2) Crecimiento lineal: Existe \(K > 0\) tal que \[ |a(t,u)| + |b(t,u)| + \left( \int_E |c(t, u, v)|^2 \phi(dv) \right)^{1/2} \leq K(1 + |u|), \tag{H2} \] para todo \(t \in [0,T]\) y \(u \in \mathbb{R}\) (Platen y Bruti-Liberati (2010), Ec. (1.9.3)).

  • (H3) Momentos finitos de saltos: \(\int_E |c(t, u, v)|^4 \phi(dv) < \infty\) para todo \(t \in [0,T]\) y \(u \in \mathbb{R}\) (Wang et al. (2007), Hipótesis 3.1).

  • (H4) Estabilidad numérica: Existe \(\Delta t_0 > 0\) y \(\delta > 0\) tales que \(1 - r(t_{n+1})\Delta t \geq \delta > 0\) para todo \(\Delta t \leq \Delta t_0\) (Higham y Kloeden (2005), Lema 2.1).

  • (H5)Condiciones iniciales: \(\mathbb{E}[|U_0|^2] < \infty\) y \(U_0 = _0\).

Bajo las hipótesis (H1)–(H5), existe una constante \(C > 0\), independiente de \(\Delta t\), tal que: \[\max_{0 \leq n \leq N} \left( \mathbb{E}\left[ |U_{t_n} - U_n|^2 \right] \right)^{1/2} \leq C \sqrt{\Delta t}. \tag{6.39}\]

Demostración. Reescribimos la SDEJ(Ecuación 6.2) en forma estándar (Platen y Bruti-Liberati (2010), Ec. (1.8.2)):

\[dU_t = a(t,U_t)\,dt + b(t,U_t)\,dW_t + \int_E c(t,U_{t-},v)\,p_\phi(dv, dt), \tag{6.40}\]

donde \(a(t,u) = c(t) + r(t)u\) y \(b(t,u) = \sigma(t)u\). Las hipótesis (H1) y (H2) satisfacen las condiciones del de Platen y Bruti-Liberati (2010), que garantiza la existencia y unicidad de una solución fuerte \(U_t\). Además asegura que la solución tiene momentos finitos de segundo orden: \[\mathbb{E}\left[ \sup_{0 \leq t \leq T} |U_t|^2 \right] \leq C_1 (1 + \mathbb{E}[|U_0|^2]) < \infty. \tag{6.41}\] Esta cota es esencial para controlar el crecimiento de la solución exacta durante la propagación del error.

Definimos el error local \(R_n = U_{t_{n+1}} - U_{n+1}\) y lo descomponemos como: \[R_n = R_n^{\text{dif}} + R_n^{\text{jump}} + R_n^{\text{implicit}}, \tag{6.42}\]

donde \(R_n^{\text{dif}}\) es el error de difusión, \(R_n^{\text{jump}}\) el error de saltos y \(R_n^{\text{implicit}}\) el error de implicitidad.

Por la expansión de Euler-Maruyama para difusiones puras (Higham y Kloeden (2005), Ec. (2.3)), el error local de la parte difusiva satisface: \[\mathbb{E}\left[ |R_n^{\text{dif}}|^2 \mid \mathcal{F_{t_n}} \right] \leq K_1 (1 + |U_n|^2) (\Delta t)^2 \tag{6.43}\]

Esta cota es de orden \((\Delta t)^3\) en \(L^2\) debido a que el esquema semi-implícito utiliza una aproximación de Euler para la parte difusiva, cuyo error local en \(L^2\) es \(\mathcal{O}((\Delta t)^{3/2})\). El error local asociado a los saltos es: \[\begin{aligned} R_n^{\text{jump}} = {}& \int_{t_n}^{t_{n+1}} \int_E c(s, U_{s-}, v)\,\tilde{p}_\phi(dv, ds) \\ & - \sum_{i=1}^{\Delta N_n} c(t_n, U_n, \xi_i) \\ & + \int_E c(t_n, U_n, v) \phi(dv)\,\Delta t. \end{aligned} \tag{6.44}\]

Aplicando la desigualdad de Cauchy-Schwarz, la condición de Lipschitz (H1) y las propiedades de la medida compensada (Platen y Bruti-Liberati (2010), Ec. (4.2.1)), obtenemos:

\[\begin{align} \mathbb{E}\left[ |R_n^{\text{jump}}|^2 \mid \mathcal{F_{t_n}} \right] &\leq 2\mathbb{E}\left[ \left| \int_{t_n}^{t_{n+1}} \int_E [c(s, U_{s-}, v) - c(t_n, U_n, v)] \tilde{p}_\phi(dv, ds) \right|^2 \mid \mathcal{F_{t_n}} \right] \nonumber \\ &\quad \ \ \ \ \ \ \ \ \ \ \ + 2\mathbb{E}\Biggl[ \Biggl| \int_{t_n}^{t_{n+1}} \int_E c(t_n, U_n, v) \tilde{p}_\phi(dv, ds) - \sum_{i=1}^{\Delta N_n} c(t_n, U_n, \xi_i) \\ &\quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ + \int_E c(t_n, U_n, v) \phi(dv)\Delta t \Biggl|^2 \mid \mathcal{F_{t_n}} \Biggr] \nonumber \\ &\leq 2L_c^2 \mathbb{E}\left[ \int_{t_n}^{t_{n+1}} |U_s - U_n|^2 ds \mid \mathcal{F_{t_n}} \right] \int_E \phi(dv) + 2 \int_E |c(t_n, U_n, v)|^2 \phi(dv) \Delta t \nonumber \\ &\leq K_2 \Delta t + K_3 |U_{t_n} - U_n|^2 \Delta t, \end{align} \tag{6.45}\] donde en la última desigualdad usamos la cota de momentos de la solución exacta Ecuación 6.41 y la condición de crecimiento lineal (H2). El primer término proviene de la variación del estado durante el intervalo, y el segundo de la diferencia entre el estado exacto y el numérico. Expandiendo el denominador mediante serie de Taylor (Wang et al. (2007), Lema 3.1): \[\frac{1}{1 - r(t_{n+1})\Delta t} = 1 + r(t_{n+1})\Delta t + \mathcal{O}((\Delta t)^2), \tag{6.46}\]

lo que implica que el error introducido por la implicitidad satisface: \[\mathbb{E}\left[|R_n^{\text{implicit}}|^2 \mid \mathcal{F_{t_n}} \right] \leq K_4 (1 + |U_n|^2) (\Delta t)^4. \tag{6.47}\]

Este término es de orden superior y no afectará el orden global de convergencia.

Combinando Ecuación 6.43, Ecuación 6.45 y Ecuación 6.47, y usando \((a+b+c)^2 \leq 3a^2 + 3b^2 + 3c^2\): \[\mathbb{E}\left[ |R_n|^2 \mid \mathcal{F_{t_n}} \right] \leq K_5 (1 + |U_n|^2 + |e_n|^2) \Delta t, \tag{6.48}\] donde \(e_n = U_{t_n} - U_n\) es el error global y usamos que \(\Delta t \leq 1\) para absorber los términos de orden superior. Definiendo \(M_n = \mathbb{E}\left[ \sup_{0 \leq k \leq n} |U_k|^2 \right]\), y tomando cuadrados en Ecuación 6.38 seguido de esperanza, aplicamos las cotas de Lipschitz (H1), crecimiento lineal (H2) y la cota de estabilidad (H4). Siguiendo el argumento del Teorema 3.2 de Higham y Kloeden (2005), obtenemos: \[M_{n+1} \leq (1 + K_6 \Delta t) M_n + K_7 \Delta t. \tag{6.49}\] Iterando esta desigualdad y usando \(M_0 = \mathbb{E}[|U_0|^2] < \infty\): \[M_N \leq M_0 e^{K_6 T} + \frac{K_7}{K_6} (e^{K_6 T} - 1) < \infty. \tag{6.50}\] Por lo tanto, el esquema numérico tiene momentos finitos uniformes, condición esencial para la convergencia fuerte. Para la descomposición del error global definimos \(e_n = U_{t_n} - U_n\). Usando la fórmula de Itô para la solución exacta y la definición del esquema Ecuación 6.38, obtenemos: \[e_{n+1} = e_n + \delta a_n \Delta t + \delta b_n \Delta W_n + \delta c_n + R_n, \tag{6.51}\] donde: \[\begin{aligned} \delta a_n &= a(t_n, U_{t_n}) - a(t_n, U_n), \\ \delta b_n &= b(t_n, U_{t_n}) - b(t_n, U_n), \\ \delta c_n &= \int_E [c(t_n, U_{t_n}, v) - c(t_n, U_n, v)] \phi(dv) \Delta t. \end{aligned}\] Por la condición de Lipschitz (H1): \[|\delta a_n| \leq L_a |e_n|, \quad |\delta b_n| \leq L_b |e_n|, \quad |\delta c_n| \leq L_c |e_n| \Delta t. \tag{6.52}\] Tomamos cuadrados en Ecuación 6.51 y esperanza condicional \(\mathbb{E}[\cdot \mid \mathcal{F_{t_n}}].\) Usando \(\mathbb{E}[\Delta W_n \mid \mathcal{F_{t_n}}] = 0\), \(\mathbb{E}[(\Delta W_n)^2 \mid \mathcal{F_{t_n}}] = \Delta t\), y que los términos con \(\Delta W_n\) tienen esperanza cero: \[\begin{aligned} \mathbb{E}[|e_{n+1}|^2 \mid \mathcal{F_{t_n}}] &\leq |e_n|^2 + 2e_n \cdot \delta a_n \Delta t + 2e_n \cdot \delta c_n + |\delta b_n|^2 \Delta t + \mathbb{E}[|R_n|^2 \mid \mathcal{F_{t_n}}] \nonumber \\ &\leq |e_n|^2 + 2L_a |e_n|^2 \Delta t + 2L_c |e_n|^2 \Delta t + L_b^2 |e_n|^2 \Delta t + K_5 (1 + |U_n|^2 + |e_n|^2) \Delta t \nonumber \\ &\leq (1 + K_8 \Delta t) |e_n|^2 + K_9 (1 + |U_n|^2) \Delta t, \end{aligned} \tag{6.53}\] donde \(K_8 = 2L_a + 2L_c + L_b^2 + K_5\) y \(K_9 = K_5\). Usando la desigualdad \(|U_n|^2 \leq 2|U_{t_n}|^2 + 2|e_n|^2\) y la cota de momentos de la solución exacta Ecuación 6.41: \[\mathbb{E}[|e_{n+1}|^2 \mid \mathcal{F_{t_n}}] \leq (1 + K_{10} \Delta t) |e_n|^2 + K_{11} \Delta t,\] donde \(K_{10} = K_8 + 2K_9\) y \(K_{11} = 2K_9 C_1 (1 + \mathbb{E}[|U_0|^2]) + K_9\). Tomando esperanza total y definiendo \(E_n = \mathbb{E}[|e_n|^2]\): \[E_{n+1} \leq (1 + K_{10} \Delta t) E_n + K_{11} \Delta t. \tag{6.54}\] Iteramos la desigualdad Ecuación 6.54 con \(E_0 = 0\): \[E_N \leq K_{11} \Delta t \sum_{k=0}^{N-1} (1 + K_{10} \Delta t)^k = K_{11} \Delta t \frac{(1 + K_{10} \Delta t)^N - 1}{K_{10} \Delta t}. \tag{6.55}\] Usando la desigualdad \((1 + x)^N \leq e^{Nx}\) para \(x \geq 0\) y \(N\Delta t = T\) Platen y Bruti-Liberati (2010), (Lema 1.2.1): \[(1 + K_{10} \Delta t)^N \leq e^{K_{10} T}. \tag{6.56}\] Sustituyendo en Ecuación 6.55: \[E_N \leq K_{11} \frac{e^{K_{10} T} - 1}{K_{10}} \Delta t = C \Delta t, \tag{6.57}\] donde \(C = K_{11} (e^{K_{10} T} - 1)/K_{10}\) es independiente de \(\Delta t\). Tomando raíz cuadrada: \[\left( \mathbb{E}[|U_T - U_N|^2] \right)^{1/2} \leq C \sqrt{\Delta t}. \quad \square\]

6.5 Método Monte Carlo

El método Monte Carlo (MC) constituye una clase de algoritmos computacionales que se basan en el muestreo aleatorio repetido para obtener resultados numéricos. En el contexto de esta tesis, el objetivo fundamental es aproximar cantidades de interés deterministas, tales como integrales multidimensionales o valores esperados de sistemas estocásticos complejos, mediante la simulación de un sistema un gran número de veces bajo condiciones aleatorias.

6.5.1 Estimador de Monte Carlo

La base teórica del método reside en la construcción de un estimador no sesgado basado en muestras independientes.

Definición 6.1 (Estimador de Monte Carlo).Sean \(X_1, X_2, \dots, X_N\) variables aleatorias independientes e idénticamente distribuidas (i.i.d.) con la misma distribución que \(X\). Definimos el estimador de Monte Carlo \(\hat{I}_N\) como la media muestral de las realizaciones de la función \(f\): \[ \hat{I}_N = \frac{1}{N} \sum_{i=1}^{N} f(X_i). \tag{6.58}\]

La validez de este estimador se sustenta en dos propiedades estadísticas fundamentales: su insesgamiento y su consistencia.

El estimador \(\hat{I}_N\) es un estimador insesgado de \(I\). Es decir, \(\mathbb{E}[\hat{I}_N] = I\).

Demostración. Dado que la esperanza es un operador lineal y las variables \(X_i\) son i.i.d.: \[\begin{aligned} \mathbb{E}[\hat{I}_N] &= \mathbb{E}\left[ \frac{1}{N} \sum_{i=1}^{N} f(X_i) \right] \\ &= \frac{1}{N} \sum_{i=1}^{N} \mathbb{E}[f(X_i)] \\ &= \frac{1}{N} \sum_{i=1}^{N} I = \frac{1}{N} \cdot N \cdot I = I. \end{aligned} \] lo que prueba el resultado. \(\quad \square\)

6.5.2 Convergencia y rigor asintótico.

Para justificar teóricamente por qué simular el sistema “miles o millones de veces” (es decir, hacer tender \(N \to \infty\)) nos permite aproximar el valor real del sistema, debemos recurrir a los teoremas límite de la teoría de la probabilidad.

6.5.2.1 Ley Fuerte de los Grandes Números.

La convergencia casi segura del estimador está garantizada por la Ley Fuerte de los Grandes Números.

Teorema 6.3 (Ley Fuerte de los Grandes Números). Sea \(\{X_i\}_{i \geq 1}\) una sucesión de variables aleatorias i.i.d. con \(\mathbb{E}[|f(X_1)|] < \infty\). Entonces, el estimador de Monte Carlo converge casi seguramente (c.s.) al valor verdadero \(I\): \[\mathbb{P}\left( \lim_{N \to \infty} \hat{I}_N = I \right) = 1.\]

Este teorema proporciona la garantía teórica de que, al aumentar el número de simulaciones, el error de aproximación tiende a cero con probabilidad unitaria. Sin embargo, la SLLN no cuantifica la tasa de convergencia ni la distribución del error para un \(N\) finito.

Para evaluar la precisión del método y construir intervalos de confianza, utilizamos el Teorema del Límite Central. Asumamos que la varianza de la función es finita, es decir, \(\sigma^2 = \mathbb{V}[f(X)] < \infty\).

Teorema 6.4 (Teorema del Límite Central para MC) Bajo las condiciones de varianza finita, la distribución del error normalizado converge en distribución a una normal estándar: \[ \sqrt{N} \frac{\hat{I}_N - I}{\sigma} \xrightarrow{d} \mathbb{N}(0, 1) \quad \text{cuando } N \to \infty.\]

De este resultado se desprende una conclusión crítica para la simulación computacional: el error estándar del estimador decae proporcionalmente a \(N^{-1/2}\), independientemente de la dimensión del problema. \[\text{Error}(\hat{I}_N) \approx \frac{\sigma}{\sqrt{N}}. \tag{6.59}\]

Esta propiedad es distintiva del método Monte Carlo. Mientras que los métodos deterministas sufren una degradación exponencial en su precisión a medida que aumenta la dimensionalidad del espacio de integración, la tasa de convergencia \(O(N^{-1/2})\) de Monte Carlo permanece invariante ante la dimensionalidad, lo que lo hace superior para sistemas complejos de alta dimensión.

6.5.3 Generación de variables aleatorias.

La implementación práctica de la ecuación Ecuación 6.58 requiere la capacidad de generar muestras \(X_i\) de una distribución de probabilidad específica \(P\). Matemáticamente, esto se logra a menudo mediante el Método de la Transformada Inversa.

Sea \(U \sim \text{Uniforme}(0,1)\) una variable aleatoria uniforme. Si \(F_X(x) = \mathbb{P}(X \leq x)\) es la función de distribución acumulada (CDF) de \(X\), y \(F_X^{-1}\) es su inversa generalizada definida como: \[F_X^{-1}(u) = \inf \{ x \in \mathbb{R} : F_X(x) \geq u \},\] entonces la variable aleatoria \(Y = F_X^{-1}(U)\) tiene la misma distribución que \(X\).

En sistemas complejos donde \(F_X^{-1}\) no tiene forma cerrada, se emplean métodos de rechazo (Rejection Sampling) o cadenas de Markov (MCMC), los cuales garantizan que la distribución estacionaria de la cadena generada coincida con la distribución objetivo \(P\).

6.5.4 Técnicas de Reducción de Varianza.

Dado que la convergencia es lenta (\(1/\sqrt{N}\)), reducir la varianza \(\sigma^2\) es tan importante como aumentar \(N\). Una varianza menor implica que se requieren menos simulaciones para alcanzar un nivel de precisión deseado.

6.5.4.1 Variables de control.

Una de las técnicas más rigurosas es el uso de variables de control. Supongamos que existe una variable aleatoria \(Y\) correlacionada con \(f(X)\) cuya esperanza \(\mu_Y = \mathbb{E}[Y]\) es conocida analíticamente. Definimos un nuevo estimador: \[\hat{I}_{CV}(c) = \hat{I}_N - c(\hat{Y}_N - \mu_Y),\] donde \(\hat{Y}_N = \frac{1}{N}\sum Y_i\) y \(c\) es una constante. Este estimador sigue siendo insesgado para cualquier \(c \in \mathbb{R}\). La varianza de este nuevo estimador es: \[\mathbb{V}[\hat{I}_{CV}(c)] = \mathbb{V}[\hat{I}_N] + c^2 \mathbb{V}[\hat{Y}_N] - 2c \text{Cov}(\hat{I}_N, \hat{Y}_N).\]

Minimizando esta expresión respecto a \(c\), obtenemos el coeficiente óptimo: \[c^* = \frac{\text{Cov}(f(X), Y)}{\mathbb{V}(Y)}.\] Sustituyendo \(c^*\), se demuestra que la varianza se reduce en un factor de \((1 - \rho^2)\), donde \(\rho\) es el coeficiente de correlación entre \(f(X)\) e \(Y\). Si existe una alta correlación, la eficiencia del método mejora drásticamente.