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:

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

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.3}\]

\(\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:

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.\]
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}\).
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}\).
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.4}\] 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.5}\]

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.6}\]

El esquema Ecuación 6.3 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.

El error local \(R_n = X_{t_{n+1}} - Y_{n+1}\) puede descomponerse como:

\[R_n = R_n^{\text{dif}} + R_n^{\text{jump}}, \tag{6.7}\] donde

\(R_n^{\text{dif}}\) es el error de la parte difusiva,
\(R_n^{\text{jump}}\) es el error de la parte de saltos.

Por la expansión de Wagner-Platen para difusiones puras (Platen y Bruti-Liberati (2010), Ec. (5.3.19)), 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 + |Y_n|^2) (\,dt)^3. \tag{6.8}\] Para la parte de saltos, el error local es: \[R_n^{\text{jump}} = \int_{t_n}^{t_{n+1}} \int_E R(s, X_{s-}, v)\,\tilde{p}_\phi(dv, ds) - \sum_{i=1}^{\Delta N_n} R(t_n, Y_n, \xi_i) - \int_E R(t_n, Y_n, v)\phi(dv)\,dt \tag{6.9}\]

Usando propiedades de la medida de Poisson compensada (ver Platen y Bruti-Liberati (2010), Ec. (4.2.1)) y la condición de Lipschitz para \(c\), se tiene: \[\mathbb{E}\left[ |R_n^{\text{jump}}|^2 \mid \mathcal{F}_{t_n} \right] \leq K_2 \,dt + K_3 |X_{t_n} - Y_n|^2 \,dt. \tag{6.10}\] Combinando Ecuación 6.8 y Ecuación 6.10 usando la desigualdad \((a+b)^2 \leq 2a^2 + 2b^2\), obtenemos: \[ \begin{aligned} \mathbb{E}\left[ |R_n|^2 \mid \mathcal{F}_{t_n} \right] \leq {}& 2K_1 (1 + |Y_n|^2) (\,dt)^3 + 2K_2 \,dt + 2K_3 |e_n|^2 \,dt \\ & \leq K_4 (1 + |Y_n|^2 + |e_n|^2) (\,dt)^2, \end{aligned} \tag{6.11}\]

donde en el último paso usamos que \(\, dt \leq 1\) (paso de tiempo pequeño) y \(e_n = X_{t_n} - Y_n\) es el error global.

Para demostrar que el esquema tiene momentos finitos, usamos un argumento de tipo Grönwall discreto. Definimos: \[M_n = \mathbb{E}\left[ \sup_{0 \leq k \leq n} |Y_k|^2 \right]. \tag{6.12}\]

Tomando cuadrados en Ecuación 6.3 y esperanza, y usando las cotas de Lipschitz y crecimiento lineal, se obtiene: \[M_{n+1} \leq (1 + K_5 \, dt) M_n + K_6 \, dt. \tag{6.13}\] Iterando esta desigualdad y usando que \(M_0 = \mathbb{E}[|Y_0|^2] < \infty\), se obtiene: \[M_N \leq M_0 e^{K_5 T} + \frac{K_6}{K_5} (e^{K_5 T} - 1) < \infty. \tag{6.14}\] Por lo tanto, el esquema numérico tiene momentos finitos uniformes.

Definimos el error global \(e_n = X_{t_n} - Y_n\). Usando la fórmula de Itô para la solución exacta y la definición del esquema, podemos escribir: \[e_{n+1} = e_n + \delta a_n \,dt + \delta b_n \,dW_n + \delta c_n + R_n, \tag{6.15}\] donde

\(\delta a_n = a(t_n, X_{t_n}) - a(t_n, Y_n)\),
\(\delta b_n = b(t_n, X_{t_n}) - b(t_n, Y_n)\),
\(\delta c_n = \int_E [c(t_n, X_{t_n}, v) - c(t_n, Y_n, v)] \phi(\,D v) \, dt\),
\(R_n\) es el residuo total (error local).

Por la condición de Lipschitz: \[|\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| \, dt, \tag{6.16}\]

donde \(L_a, L_b, L_c\) son las constantes de Lipschitz de \(a\), \(b\) y \(c\), respectivamente.

Tomamos cuadrados en Ecuación 6.15: \[\begin{aligned} |e_{n+1}|^2 &= |e_n|^2 + 2e_n \cdot \delta a_n \, dt + 2e_n \cdot \delta b_n \, dW_n + 2e_n \cdot \delta c_n \\ &\quad + |\delta a_n|^2 (\, dt)^2 + 2\delta a_n \cdot \delta b_n \, dt \, dW_n + 2\delta a_n \cdot \delta c_n \, dt \\ &\quad + |\delta b_n|^2 (\, dW_n)^2 + 2\delta b_n \cdot \delta c_n \, dW_n + |\delta c_n|^2 \\ &\quad + 2e_n \cdot R_n + 2\delta a_n \cdot R_n \, dt + 2\delta b_n \cdot R_n \, dW_n + 2\delta c_n \cdot R_n + |R_n|^2. \end{aligned} \tag{6.17}\] Tomamos esperanza condicional dado \(\mathcal{F}_{t_n}\). Notamos que:

\(\mathbb{E}[\,dW_n \mid \mathcal{F}_{t_n}] = 0\),
\(\mathbb{E}[(\,dW_n)^2 \mid \mathcal{F}_{t_n}] = dt\),
Los términos con \(\,dW_n\) tienen esperanza cero.

Por lo tanto: \[\begin{aligned} \mathbb{E}[|e_{n+1}|^2 \mid \mathcal{F}_{t_n}] &\leq |e_n|^2 + 2e_n \cdot \delta a_n \, dt + 2e_n \cdot \delta c_n + |\delta b_n|^2 \, dt \\ &\quad + |\delta a_n|^2 (\, dt)^2 + 2|\delta a_n| |\delta c_n| \, dt + |\delta c_n|^2 \\ &\quad + 2|e_n| |R_n| + 2|\delta a_n| |R_n| \, dt + 2|\delta c_n| |R_n| + |R_n|^2. \end{aligned} \tag{6.18}\]

Usando las cotas de Lipschitz Ecuación 6.16 y la desigualdad \(2ab \leq a^2 + b^2\), obtenemos: \[\begin{aligned} \mathbb{E}[|e_{n+1}|^2 \mid \mathcal{F}_{t_n}] &\leq |e_n|^2 + 2L_a |e_n|^2 \, dt + 2L_c |e_n|^2 \, dt + L_b^2 |e_n|^2 \, dt \\ &\quad + L_a^2 |e_n|^2 (\, dt)^2 + 2L_a L_c |e_n|^2 (\, dt)^2 + L_c^2 |e_n|^2 (\, dt)^2 \\ &\quad + |e_n|^2 + |R_n|^2 + L_a^2 |e_n|^2 (\, dt)^2 + |R_n|^2 (\, dt)^2 \\ &\quad + L_c^2 |e_n|^2 (\, dt)^2 + |R_n|^2 + |R_n|^2. \end{aligned} \tag{6.19}\]

Simplificando y usando que \(\, dt \leq 1\): \[\mathbb{E}[|e_{n+1}|^2 \mid \mathcal{F}_{t_n}] \leq (1 + K_7 \, dt) |e_n|^2 + K_8 |R_n|^2, \tag{6.20}\] donde \(K_7 = 2L_a + 2L_c + L_b^2 + L_a^2 + 2L_a L_c + L_c^2 + 1\) y \(K_8 = 4 + (\, dt)^2 \leq 5\). Usando la cota del error local Ecuación 6.11: \[\mathbb{E}[|e_{n+1}|^2 \mid \mathcal{F}_{t_n}] \leq (1 + K_7 \, dt) |e_n|^2 + K_8 K_4 (1 + |Y_n|^2 + |e_n|^2) (\, dt)^2.\] Usando la desigualdad \(|Y_n|^2 \leq 2|X_{t_n}|^2 + 2|e_n|^2\) y tomando esperanza total: \[\mathbb{E}[|e_{n+1}|^2] \leq (1 + K_7 \, dt + 2K_8 K_4 (\, dt)^2) \mathbb{E}[|e_n|^2] + 2K_8 K_4 (\, dt)^2 \mathbb{E}[|X_{t_n}|^2] + K_8 K_4 (\, dt)^2. \tag{6.21}\] Definiendo \(E_n = \mathbb{E}[|e_n|^2]\) y usando la cota de momentos de la solución exacta Ecuación 6.6:

\[E_{n+1} \leq (1 + K_9 \, dt) E_n + K_{10} (\, dt)^2, \tag{6.22}\]

donde \(K_9 = K_7 + 2K_8 K_4 \, dt \leq K_7 + 2K_8 K_4 \ \text{y} K_{10} = 2K_8 K_4 C_1 (1 + \mathbb{E}[|X_0|^2]) + K_8 K_4\).

Iteramos la desigualdad Ecuación 6.22: \[\begin{aligned} E_1 &\leq (1 + K_9 \, dt) E_0 + K_{10} (\, dt)^2, \\ E_2 &\leq (1 + K_9 \, dt)^2 E_0 + K_{10} (\, dt)^2 \left[1 + (1 + K_9 \, dt)\right], \\ &\vdots \\ E_N &\leq (1 + K_9 \, dt)^N E_0 + K_{10} (\, dt)^2 \sum_{k=0}^{N-1} (1 + K_9 \, dt)^k. \end{aligned} \tag{6.23}\]

Como \(E_0 = 0\) (condiciones iniciales iguales), el primer término se anula. La suma geométrica se puede acotar: \[\sum_{k=0}^{N-1} (1 + K_9 \, dt)^k = \frac{(1 + K_9 \, dt)^N - 1}{K_9 \, dt}. \tag{6.24}\]

Usando la desigualdad \((1 + x)^N \leq e^{Nx}\) para \(x \geq 0\), y notando que \(N\, dt = T\), obtenemos: \[(1 + K_9 \, dt)^N \leq e^{K_9 T} \tag{6.25}\] donde \(C = K_{10} (e^{K_9 T} - 1)/K_9\) es una constante independiente de \(\,dt\). Tomando raíz cuadrada, obtenemos finalmente: \[\left( \mathbb{E}[|X_T - Y_N|^2] \right)^{1/2} \leq C \sqrt{\,dt}. \quad \square \tag{6.26}\]

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.27}\]

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.28}\] 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.29}\]

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.30}\]

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.31}\] 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.32}\]

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.33}\]

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.34}\]

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.35}\] donde en la última desigualdad usamos la cota de momentos de la solución exacta Ecuación 6.31 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.36}\] 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.\]

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

Combinando Ecuación 6.33, Ecuación 6.35 y ?eq-implicit_error_bound_semiimplicit, 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.37}\] 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.28 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.38}\] 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.39}\] 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.28, obtenemos: \[e_{n+1} = e_n + \delta a_n \Delta t + \delta b_n \Delta W_n + \delta c_n + R_n, \tag{6.40}\] 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.41}\] Tomamos cuadrados en Ecuación 6.40 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.42}\] 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.31: \[\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.43}\] Iteramos la desigualdad Ecuación 6.43 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.44}\] 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.45}\] Sustituyendo en Ecuación 6.44: \[E_N \leq K_{11} \frac{e^{K_{10} T} - 1}{K_{10}} \Delta t = C \Delta t, \tag{6.46}\] 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.47}\]

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.48}\]

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.47 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.