8  Análisis de Resultados.

En este capítulo se presentan y discuten los resultados numéricos obtenidos de la simulación numérica del modelo de Cramér-Lundberg extendido con difusión y saltos, descrito en la Sección 5.2.2, propuesto para la ecuación diferencial estocástica con saltos (EDE) que modela las reservas de la aseguradora. El objetivo central es validar la robustez computacional de los métodos Milstein Compensado (Ecuación 8.1) y Semi-Implícito Compensado (Ecuación 8.2), evaluando su desempeño bajo las condiciones de mercado del sector asegurador de vivienda en la Ciudad de México.

8.1 Simulación Milstein sin Compensación.

Aunque ambos esquemas comparten el mismo orden de convergencia fuerte (\(\gamma=1/2\)) y satisfacen las condiciones de Lipschitz y crecimiento lineal, la versión compensada fue seleccionada por razones técnicas, numéricas y regulatorias alineadas con el objetivo actuarial del trabajo:

1. Menor constante de error y cancelación de sesgo:
   La compensación introduce explícitamente \(-\int_E R(t_n,Y_n,v)\phi(dv)\Delta t\) en la deriva. Esto garantiza \(\mathbb{E}[\text{término de saltos}] = 0\), eliminando el sesgo de primer orden y reduciendo la constante \(C\) en la cota de convergencia. En la práctica, el error RMS entre mallas se mantiene acotado incluso con \(\Delta t\) relativamente grandes.

2. Consistencia con la dinámica actuarial del modelo extendido:
   La Ec. (3.4) de la tesis define la deriva neta como \(\mu(t,X_t) = c(t) + r(t)X_t - \lambda_{\text{meteo}}(t)\mathbb{E}[Y^{\text{meteo}}] - \lambda_{\text{sismo}}\mathbb{E}[Y^{\text{sismo}}]\). El esquema compensado preserva esta estructura a nivel discreto, asegurando que la simulación refleje fielmente el principio de equivalencia actuarial (primas = siniestros esperados + margen).

3. Estabilidad numérica con deriva lineal \(r(t)X_t\):
   Al combinar Milstein compensado con el esquema semi-implícito (Ec. 43), el denominador \(1 - r(t_{n+1})\Delta t\) opera sobre un ruido de saltos de media cero. Esto evita amplificaciones explosivas cuando \(r(t) > 0\) (caso real en México con \(r \approx 8.9\%\) anual). El esquema no compensado arrastra la media de los saltos a la deriva, pudiendo violar la condición de estabilidad \(1 - r\Delta t > 0\) en periodos de alta siniestralidad.

4. Eficiencia en inferencia Monte Carlo y cumplimiento regulatorio:
   La CNSF exige proyecciones de reservas robustas y no sesgadas. La versión compensada reduce la varianza intrínseca del estimador, permitiendo alcanzar el IC 90% con \(M=1000\) trayectorias (como se reporta en la Sección 3.6). Además, facilita la aplicación de técnicas de reducción de varianza (variables de control) al mantener la estructura de martingala.

El código Milstein proporcionado implementado con \(a(t,x)=\mu x\), \(b(t,x)=\sigma x\) y saltos exponenciales restados directamente. Al comparar con la teoría y los resultados, se identifican limitaciones críticas.

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt

# Configuración visual opcional
plt.style.use('seaborn-v0_8-whitegrid')
def milstein_with_jumps(X0, mu, sigma, lam, beta, T, N):
    """
    Simula una trayectoria de dX = mu*X dt + sigma*X dW - dJ
    usando el método de Milstein para la parte difusiva y
    saltos de Poisson compuesto.

    Parámetros:
        X0 (float): valor inicial
        mu (float): tasa de crecimiento
        sigma (float): volatilidad
        lam (float): intensidad del proceso de Poisson (λ)
        beta (float): parámetro de la distribución exponencial de los saltos (media = 1/β)
        T (float): tiempo final
        N (int): número de pasos temporales

    Retorna:
        X (np.array): trayectoria simulada (longitud N+1)
        t (np.array): vector de tiempos
        jump_times (list): lista de tiempos donde ocurrieron saltos
        jump_sizes (list): lista de tamaños de los saltos
    """
    dt = T / N
    X = np.empty(N + 1)
    X[0] = X0

    # Para registrar los saltos (opcional, para visualización)
    jump_times = []
    jump_sizes = []

    # Pre-generar todos los incrementos brownianos
    dW = np.sqrt(dt) * np.random.randn(N)

    for i in range(N):
        # 1. Simular el número de saltos en este intervalo
        n_jumps = np.random.poisson(lam * dt)

        # 2. Calcular el tamaño total del salto
        if n_jumps > 0:
            # Tamaños de salto exponenciales
            jumps = np.random.exponential(scale=1.0/beta, size=n_jumps)
            total_jump = np.sum(jumps)
            # Registrar para visualización (tiempo al final del intervalo)
            jump_times.extend([ (i+1) * dt ] * n_jumps)
            jump_sizes.extend(jumps.tolist())
        else:
            total_jump = 0.0

        # 3. Aplicar el esquema de Milstein con el salto
        X[i+1] = (X[i]
                  + mu * X[i] * dt
                  + sigma * X[i] * dW[i]
                  + 0.5 * sigma**2 * X[i] * (dW[i]**2 - dt)
                  - total_jump)

        # Evitar que el capital se vuelva negativo (opcional, depende del modelo)
        # X[i+1] = max(X[i+1], 0.0)

    t = np.linspace(0, T, N + 1)
    return X, t, jump_times, jump_sizes
# -----------------------------
# Parámetros del modelo
# -----------------------------
T = 1.0          # Tiempo final
X0 = 100.0       # Capital inicial
mu = 0.5         # Tasa de crecimiento
sigma = 0.3      # Volatilidad
lam = 5.0        # Intensidad de saltos (λ)
beta = 0.1       # Parámetro exponencial (media del salto = 10)

# Mallas
N_gruesa = 256   # baja frecuencia
N_fina   = 1024  # alta frecuencia

# Fijar semilla para reproducibilidad
np.random.seed(123)

# Simulaciones
X_g, t_g, jt_g, js_g = milstein_with_jumps(X0, mu, sigma, lam, beta, T, N_gruesa)
X_f, t_f, jt_f, js_f = milstein_with_jumps(X0, mu, sigma, lam, beta, T, N_fina)

# Interpolar la trayectoria fina en los puntos de la gruesa
X_f_interp = np.interp(t_g, t_f, X_f)

# Calcular error absoluto
error_abs = np.abs(X_g - X_f_interp)
# -----------------------------
# Gráficos
# -----------------------------
plt.figure(figsize=(14, 10))

# Trayectorias
plt.subplot(2, 2, 1)
plt.plot(t_g, X_g, 'o-', markersize=3, label=f'Malla gruesa (N={N_gruesa})', alpha=0.8)
plt.plot(t_f, X_f, '-', linewidth=1, label=f'Malla fina (N={N_fina})', alpha=0.9)
plt.title('Milstein con Saltos: Comparación de mallas')
plt.xlabel('Tiempo $t$')
plt.ylabel('$X(t)$')
plt.legend()
plt.grid(True)

# Error absoluto
plt.subplot(2, 2, 2)
plt.plot(t_g, error_abs, 'r.-', markersize=4)
plt.title('Error absoluto (gruesa vs fina interpolada)')
plt.xlabel('Tiempo $t$')
plt.ylabel('|Error|')
plt.yscale('log')
plt.grid(True)

# Zoom en la malla fina con saltos marcados
plt.subplot(2, 1, 2)
plt.plot(t_f, X_f, '-', linewidth=1.2, label='Trayectoria fina')
# Marcar los saltos como líneas verticales
for jt in jt_f:
    plt.axvline(x=jt, color='red', linestyle='--', alpha=0.3, linewidth=0.7)
plt.title('Trayectoria fina con saltos marcados (líneas verticales rojas)')
plt.xlabel('Tiempo $t$')
plt.ylabel('$X(t)$')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

# Imprimir métricas
print(f"Error máximo absoluto: {error_abs.max():.6e}")
print(f"Error RMS:             {np.sqrt(np.mean(error_abs**2)):.6e}")
print(f"Número de saltos (gruesa): {len(jt_g)}")
print(f"Número de saltos (fina):   {len(jt_f)}")

Error máximo absoluto: 4.632973e+01
Error RMS:             1.916853e+01
Número de saltos (gruesa): 5
Número de saltos (fina):   3

La gráfica muestra la evolución de las reservas de una aseguradora de vivienda en la CDMX simulada mediante el método de Milstein con saltos compensados durante 2 años (2024-2025).

import numpy as np
import matplotlib.pyplot as plt

# -----------------------------
# Parámetros iniciales
# -----------------------------
T_total = 2  # años
N_months = 24
days_per_month = 30
N_days = N_months * days_per_month
dt = 1 / 365  # paso diario (en años)

# -----------------------------
# Funciones temporales (igual que antes)
# -----------------------------
c_monthly = np.array([
    162, 162, 163, 164, 165, 166,
    167, 168, 169, 170, 171, 172,
    173, 173, 174, 174, 175, 175,
    176, 176, 177, 177, 178, 178
])
c_daily = np.repeat(c_monthly / 30, days_per_month)

r_ILS = np.array([
    0.75, 0.68, 0.72, 0.65, 0.80, 0.95,
    0.60, 1.20, 0.50, 0.70, 0.65, 0.85,
    0.78, 0.70, 0.75, 0.82, 0.90, 1.10,
    0.65, 0.55, 0.95, 0.72, 0.68, 0.90
]) / 100

r_base = 0.0065
L = np.maximum(r_ILS - r_base, 0)
r_daily = np.repeat((1 + r_ILS) ** (1/30) - 1, days_per_month)
sigma_base = 0.02
alpha = 5
sigma_daily = sigma_base * (1 + alpha * np.repeat(L, days_per_month))

seasonal_factor = np.array([
    0.6, 0.5, 0.6, 0.7, 0.9, 1.3,
    1.4, 1.8, 1.7, 1.5, 1.0, 0.7
])
seasonal_factor_24 = np.tile(seasonal_factor, 2)
lambda_annual = 12
lambda_monthly = lambda_annual * seasonal_factor_24 / 12
lambda_daily = np.repeat(lambda_monthly, days_per_month) / 30

K = 14.5
sigma_Y = 0.5
mu_Y = np.log(np.maximum(K * np.repeat(L, days_per_month) / lambda_daily, 1e-6)) - sigma_Y**2 / 2

# -----------------------------
# Simulación Monte Carlo
# -----------------------------
np.random.seed(42)
N_sim = 1000
U_all = np.zeros((N_sim, N_days + 1))
U_all[:, 0] = 10000.0  # Reserva inicial

for sim in range(N_sim):
    for i in range(N_days):
        dW = np.sqrt(dt) * np.random.randn()
        n_jumps = np.random.poisson(lambda_daily[i])
        total_jump = np.sum(np.random.lognormal(mean=mu_Y[i], sigma=sigma_Y, size=n_jumps)) if n_jumps > 0 else 0.0

        # Método de Milstein completo
        drift = c_daily[i] * dt + r_daily[i] * U_all[sim, i] * dt
        diffusion = sigma_daily[i] * U_all[sim, i] * dW
        milstein_correction = 0.5 * (sigma_daily[i]**2) * U_all[sim, i] * (dW**2 - dt)

        U_all[sim, i+1] = U_all[sim, i] + drift + diffusion + milstein_correction - total_jump

# Calcular estadísticas
median_path = np.median(U_all, axis=0)
p05 = np.percentile(U_all, 5, axis=0)
p95 = np.percentile(U_all, 95, axis=0)

# Probabilidad de ruina
ruin_occurred = (U_all <= 0).any(axis=1)
prob_ruin_final = np.mean(U_all[:, -1] <= 0)
prob_ruin_anytime = np.mean(ruin_occurred)

print(f"Probabilidad de ruina final: {prob_ruin_final:.2%}")
print(f"Probabilidad de ruina en algún momento: {prob_ruin_anytime:.2%}")

# -----------------------------
# Visualización mejorada
# -----------------------------
time_days_U = np.arange(N_days + 1)
months_labels = [f"{m//12+2024}-{(m%12)+1:02d}" for m in range(N_months)]
month_ticks = np.arange(0, N_days, 30)

plt.figure(figsize=(14, 8))

# Gráfico principal
plt.plot(time_days_U, median_path, color='navy', linewidth=2.5, label='Mediana')
plt.fill_between(time_days_U, p05, p95, color='steelblue', alpha=0.3, label='Percentiles 5%-95%')

# Graficar 50 trayectorias de muestra en un color claro
sample_paths = 50
for i in range(sample_paths):
    plt.plot(time_days_U, U_all[i, :], color='lightsteelblue', alpha=0.6, linewidth=0.8)

plt.title('Reserva de la aseguradora — Simulación Monte Carlo (1,000 trayectorias)', fontsize=14)
plt.ylabel('Reserva (USD)', fontsize=12)
plt.xlabel('Tiempo', fontsize=12)
plt.xticks(month_ticks, months_labels, rotation=45)
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

Probabilidad de ruina final: 0.00%
Probabilidad de ruina en algún momento: 0.00%

El esquema de Milstein sin compensación es teóricamente válido y converge con orden \(1/2\), pero introduce un sesgo de deriva, infla la varianza y degrada la estabilidad en presencia de saltos frecuentes o de magnitud significativa, como los modelados para siniestros hidrometeorológicos y sísmicos en la CDMX. La versión compensada conserva el mismo orden de convergencia, pero elimina el sesgo de primer orden, preserva la estructura de martingala necesaria para las cotas de error, y se alinea con la dinámica actuarial del modelo de Cramér-Lundberg extendido.

8.2 Milstein compensado.

El esquema de Milstein para la componente difusiva, combinado con compensación de saltos, se define como: \[\begin{aligned} U_{n+1}^{\text{Mil}} = U_n &+ \mu(t_n, U_n)\Delta t + \sigma(t_n, U_n)\Delta W_n \\ & + \frac{1}{2}\sigma(t_n, U_n)\sigma_U(t_n, U_n)\left[(\Delta W_n)^2 - \Delta t\right] - \Delta J_n + \mathcal{C}_n\Delta t, \end{aligned} \tag{8.1}\]

donde \(\mathcal{C}_n = -\left[\lambda_{\text{meteo},n}\mathbb{E}[Y^{\text{meteo}}_n] + \lambda_{\text{sismo},n}\mathbb{E}[Y^{\text{sismo}}_n]\right]\) es el término de compensación. Bajo las condiciones de Lipschitz y crecimiento lineal verificadaseste esquema garantiza convergencia fuerte de orden \(\mathcal{O}(\Delta t)\) para la parte continua Kloeden y Platen (1992). El término de corrección \(\frac{1}{2}\sigma \sigma_U[(\Delta W)\^2 - \Delta t]\) es esencial para capturar la no linealidad de la difusión \(\sigma(t)U(t)\), especialmente relevante ante la volatilidad estacional observada en \(\sigma_t\) (?tbl-volatilidad_total).

8.3 Semi-Implícito compensado.

El esquema semi-implícito trata el término de interés \(r(t)U(t)\) de forma implícita: \[U_{n+1}^{\text{SI}} = \frac{U_n\left[1 + \sigma(t_n, U_n)\Delta W_n\right] + c(t_n)\Delta t + \mathcal{C}_n\Delta t - \Delta J_n^{\text{bruto}}}{1 - r(t_n)\Delta t}, \tag{8.2}\]

donde \(\Delta J_n^{\text{bruto}}\) es la suma no compensada de saltos ocurridos. Este esquema es incondicionalmente estable para \(r(t) > 0\) Higham et al. (2002), ya que el denominador \(1 - r(t_n)\Delta t > 0\) (dado \(r_{\max} \approx 0.6\%\) mensual y \(\Delta t = 1/365\)). La estabilidad numérica es crítica ante la alta sensibilidad de las reservas a fluctuaciones en \(r(t)\), como se observa en la estructura de tasas de Banxico (Banco de México (2024),).

8.4 Análisis Cuantitativo de la Simulación

Se presentan las definiciones operativas de los esquemas numéricos, destacando sus fundamentos teóricos: la corrección de Milstein para la parte difusiva y el tratamiento implícito del término de interés para garantizar estabilidad incondicional.

En el cual se ejecutaron \(N_{\text{sim}} = 1000\) trayectorias Monte Carlo para el periodo histórico (enero 2024 - diciembre 2025) y pronóstico de enero 2026, con semilla fija (\(2024-2025\) para histórico, \(2026\) para pronóstico) para reproducibilidad. La cartera simulada corresponde a \(N_{\text{pólizas}} = 6,800\) pólizas, representando aproximadamente el \(40\%\) del mercado total de seguros de vivienda en la CDMX

Comparación cuantitativa de resultados para enero 2026 Milstein vs. Semi-Implícito (millones USD). {#compara_mil_semi}

Métrica	Milstein Compensado	Semi-Implícito Compensado	Diferencia Relativa
Reserva inicial (\(U_0\))	\(520.00\)	\(520.00\)	\(0.00\%\)
Mediana final	\(528.45\)	\(540.60\)	\(+2.30\%\)
Límite inferior IC \(90\%\) (\(5\%\))	\(512.30\)	\(521.50\)	\(+1.80\%\)
Límite superior IC \(90\%\) (\(95\%\))	\(545.80\)	\(559.20\)	\(+2.46\%\)
Ancho IC \(90\%\)	\(33.50\)	\(37.70\)	\(+12.54\%\)
Desviación estándar final	\(8.21\)	\(9.15\)	\(+11.45\%\)
Error estándar de la mediana	\(0.26\)	\(0.29\)	\(+11.54\%\)

import numpy as np
import plotly.graph_objects as go
# -----------------------------
# PARTE 0: CONFIGURACIÓN INICIAL
# -----------------------------
np.random.seed(2026)
tipo_cambio = 18.5
dt = 1 / 365  # Paso diario

# -----------------------------
# PARTE 1: CALIBRACIÓN DE PARÁMETROS
# -----------------------------

N_polizas = 6800
prima_base_por_poliza = 850.00
delta = 0.0015
X0_hist = 500e6  # Capital inicial histórico
U0_forecast = 520e6  # Capital inicial para pronóstico (ejemplo)

# Inflación y primas
inflacion_mensual = np.array([
    0.52, 0.38, 0.38, 0.41, 0.45, 0.49, 0.55, 0.61, 0.58, 0.47, 0.42, 0.39,
    0.35, 0.50, 0.37, 0.40, 0.44, 0.48, 0.53, 0.59, 0.56, 0.46, 0.41, 0.38
]) / 100.0

g = inflacion_mensual + delta
cum_g = np.cumsum(g)
primas_por_poliza_mxn = prima_base_por_poliza * np.exp(cum_g)
primas_totales_mxn = primas_por_poliza_mxn * N_polizas
c_t_mensual_usd = primas_totales_mxn / tipo_cambio
c_t_mensual_usd = np.round(c_t_mensual_usd, 2)

# --- Intensidades calibradas ---
lambda_meteo_annual = 6.0
lambda_sismo = 0.15 #0.08
seasonal_raw = [0.3, 0.3, 0.4, 0.6, 1.0, 1.8, 2.2, 2.5, 2.8, 2.3, 1.2, 0.5]
seasonal_norm = np.array(seasonal_raw) / np.mean(seasonal_raw)

# --- Parámetros de severidad (ahora dependientes del estado) ---
ceiling_per_policy = 5000.0
COEF_VARIACION_METEO = 0.1
COEF_VARIACION_SISMO = 0.2

# --- Rendimiento y volatilidad ---
r_t_mensual = np.array([
    0.606, 0.602, 0.600, 0.597, 0.593, 0.589,
    0.586, 0.582, 0.578, 0.574, 0.570, 0.566,
    0.562, 0.559, 0.555, 0.551, 0.547, 0.543,
    0.539, 0.535, 0.532, 0.528, 0.524, 0.520
]) / 100.0

sigma_t_mensual = np.array([
    0.44, 0.47, 0.50, 0.53, 0.59, 0.74,
    0.89, 1.34, 1.19, 0.98, 0.74, 0.68,
    0.65, 0.60, 0.58, 0.60, 0.65, 0.74,
    0.89, 0.82, 0.95, 0.86, 0.74, 0.68
]) / 100.0

# -----------------------------
# PARTE 2: FUNCIÓN AUXILIAR PARA LA COMPENSACIÓN
# -----------------------------

def expected_loss_per_event(deviation, pct_low, pct_high, n_pols):
    """
    Calcula la pérdida esperada por un evento dado el estado actual.
    """
    E_pct = (pct_low + pct_high) / 2
    n_aff_mean = n_pols * E_pct #Valor esperado de la distribución uniforme para el número de afectados
    # mu_per_policy = deviation  <-- This line is now redundant as deviation is used directly
    return n_aff_mean * deviation

# -----------------------------
# PARTE 3: SIMULACIÓN HISTÓRICA (Ene 2024 - Dic 2025) CON COMPENSACIÓN
# -----------------------------

T_hist = 2.0
N_hist = int(T_hist / dt)
N_sim_hist = 1000

dias_hist = np.linspace(0, T_hist, N_hist + 1)
meses_enteros = np.arange(24)

r_hist = np.interp(dias_hist[:-1], meses_enteros, r_t_mensual)
sigma_hist = np.interp(dias_hist[:-1], meses_enteros, sigma_t_mensual)
c_hist = np.interp(dias_hist[:-1], meses_enteros, c_t_mensual_usd)

X_hist = np.zeros((N_sim_hist, N_hist + 1))
X_hist[:, 0] = X0_hist

print("Simulando periodo histórico (Ene 2024 - Dic 2025) con Milstein y compensación...")

for sim in range(N_sim_hist):
    for i in range(N_hist):
        t_dia = dias_hist[i]
        mes_actual = int(t_dia) % 12

        lam_meteo_dia = (lambda_meteo_annual / 365) * seasonal_norm[mes_actual]
        lam_sismo_dia = lambda_sismo / 365

        n_meteo = np.random.poisson(lam_meteo_dia)
        n_sismo = np.random.poisson(lam_sismo_dia)

        x_current = X_hist[sim, i]
        deviation = (abs(x_current - X0_hist))/N_polizas
        if deviation > ceiling_per_policy:
            deviation = ceiling_per_policy

        total_jump = 0.0
        compensacion = 0.0

        # --- Eventos Meteorológicos ---
        if n_meteo > 0:
            pct_affected = np.random.uniform(0.005, 0.02)
            n_affected = max(1, int(N_polizas * pct_affected))
            mu_per_policy = deviation
            sigma_per_policy = COEF_VARIACION_METEO * mu_per_policy
            losses = np.random.normal(mu_per_policy, sigma_per_policy, size=n_affected)
            losses = np.clip(losses, 0.0, None)
            total_jump += np.sum(losses)

        # Compensación meteorológica
        E_loss_meteo = expected_loss_per_event(
            deviation, 0.005, 0.02, N_polizas
        )
        compensacion += (-1)*lam_meteo_dia * E_loss_meteo

        # --- Eventos Sísmicos ---
        if n_sismo > 0:
            pct_affected = np.random.uniform(0.01, 0.10)
            n_affected = max(1, int(N_polizas * pct_affected))
            mu_per_policy = deviation
            sigma_per_policy = COEF_VARIACION_SISMO * mu_per_policy
            losses = np.random.normal(mu_per_policy, sigma_per_policy, size=n_affected)
            losses = np.clip(losses, 0.0, None)
            total_jump += np.sum(losses)

        # Compensación sísmica
        E_loss_sismo = expected_loss_per_event(
            deviation, 0.01, 0.10, N_polizas
        )
        compensacion += (-1)*lam_sismo_dia * E_loss_sismo

        # --- Actualización con Milstein + Compensación ---
        dW = np.sqrt(dt) * np.random.randn()
        drift_base = c_hist[i] * dt + r_hist[i] * x_current * dt
        drift_compensado = drift_base + compensacion
        diffusion = sigma_hist[i] * x_current * dW
        milstein_corr = 0.5 * (sigma_hist[i]**2) * x_current * (dW**2 - dt)
        X_hist[sim, i+1] = x_current + drift_compensado + diffusion + milstein_corr - total_jump

U0_forecast = np.median(X_hist[:, -1])
print(f"Condición inicial para pronóstico (Ene 2026): ${U0_forecast/1e6:.2f}M")
# -----------------------------
# PARTE 4: SIMULACIÓN DEL PRONÓSTICO (Ene 2026) CON COMPENSACIÓN
# -----------------------------

forecast_days = 30
mes_pronostico = 0

# Promedios para enero (mes 0 y mes 12)
c_ene_avg = (c_t_mensual_usd[0] + c_t_mensual_usd[12]) / 2
r_ene_avg = (r_t_mensual[0] + r_t_mensual[12]) / 2
sigma_ene_avg = (sigma_t_mensual[0] + sigma_t_mensual[12]) / 2
sf_ene_avg = seasonal_norm[0]  # Enero es el mes 0

c_daily_fc = np.full(forecast_days, c_ene_avg / 30.0)
r_daily_fc = np.full(forecast_days, r_ene_avg)
sigma_daily_fc = np.full(forecast_days, sigma_ene_avg)
lambda_meteo_dia_fc = (lambda_meteo_annual / 365) * sf_ene_avg
lambda_sismo_dia_fc = lambda_sismo / 365

np.random.seed(2027)
N_sim_fc = 1000
U_forecast = np.zeros((N_sim_fc, forecast_days + 1))
U_forecast[:, 0] = U0_forecast

print("Simulando pronóstico para Enero 2026 con Milstein y compensación...")

for sim in range(N_sim_fc):
    for i in range(forecast_days):
        u_current = U_forecast[sim, i]
        deviation = (abs(u_current - U0_forecast))/N_polizas
        if deviation > ceiling_per_policy:
            deviation = ceiling_per_policy
        total_jump = 0.0
        compensacion = 0.0

        # --- Simular eventos ---
        n_meteo = np.random.poisson(lambda_meteo_dia_fc)
        n_sismo = np.random.poisson(lambda_sismo_dia_fc)

        # --- Eventos Meteorológicos ---
        if n_meteo > 0:
            pct_affected = np.random.uniform(0.005, 0.02)
            n_affected = max(1, int(N_polizas * pct_affected))
            mu_per_policy = deviation
            sigma_per_policy = COEF_VARIACION_METEO * mu_per_policy
            losses = np.random.normal(mu_per_policy, sigma_per_policy, size=n_affected)
            losses = np.clip(losses, 0.0, None)
            total_jump += np.sum(losses)

        # Compensación meteorológica
        E_loss_meteo = expected_loss_per_event(
            deviation, 0.005, 0.02, N_polizas
        )
        compensacion += (-1)*lambda_meteo_dia_fc * E_loss_meteo

        # --- Eventos Sísmicos ---
        if n_sismo > 0:
            pct_affected = np.random.uniform(0.01, 0.10)
            n_affected = max(1, int(N_polizas * pct_affected))
            mu_per_policy = deviation
            sigma_per_policy = COEF_VARIACION_SISMO * mu_per_policy
            losses = np.random.normal(mu_per_policy, sigma_per_policy, size=n_affected)
            losses = np.clip(losses, 0.0, None)
            total_jump += np.sum(losses)

        # Compensación sísmica
        E_loss_sismo = expected_loss_per_event(
            deviation, 0.01, 0.10, N_polizas
        )
        compensacion +=(-1)* lambda_sismo_dia_fc * E_loss_sismo

        # --- Actualización con Milstein + Compensación ---
        dW = np.sqrt(dt) * np.random.randn()
        drift_base = c_daily_fc[i] * dt + r_daily_fc[i] * u_current * dt
        drift_compensado = drift_base + compensacion
        diffusion = sigma_daily_fc[i] * u_current * dW
        milstein_corr = 0.5 * (sigma_daily_fc[i]**2) * u_current * (dW**2 - dt)
        U_forecast[sim, i+1] = u_current + drift_compensado + diffusion + milstein_corr - total_jump

# -----------------------------
# PARTE 5: CÁLCULO DE ESTADÍSTICAS
# -----------------------------

median_hist = np.median(X_hist, axis=0) / 1e6
p05_hist = np.percentile(X_hist, 5, axis=0) / 1e6
p95_hist = np.percentile(X_hist, 95, axis=0) / 1e6

median_fc = np.median(U_forecast, axis=0) / 1e6
p05_fc = np.percentile(U_forecast, 5, axis=0) / 1e6
p95_fc = np.percentile(U_forecast, 95, axis=0) / 1e6

# -----------------------------
# PARTE 6: DETECCIÓN DE SALTOS EN LA MEDIANA
# -----------------------------

def detect_jumps(median_series):
    """Detecta saltos en una serie temporal usando umbral de 3 desviaciones estándar."""
    median_usd = median_series * 1e6
    log_returns = np.log(median_usd[1:] / median_usd[:-1])
    threshold = 3 * np.std(log_returns)
    jump_flags = np.abs(log_returns) > threshold
    jump_indices = np.where(jump_flags)[0] + 1
    jump_times = np.arange(len(median_series))[jump_indices]
    jump_values = median_series[jump_indices]
    return jump_times, jump_values

# Detectar saltos
time_historical = np.arange(N_hist + 1)
jump_times_hist, jump_values_hist = detect_jumps(median_hist)

time_forecast = np.arange(N_hist, N_hist + forecast_days + 1)
jump_times_fc, jump_values_fc = detect_jumps(median_fc)

# -----------------------------
# PARTE 7: GRÁFICA INTERACTIVA CON PLOTLY
# -----------------------------

fig = go.Figure()

# Área de incertidumbre histórico
fig.add_trace(go.Scatter(
    x=np.concatenate([time_historical, time_historical[::-1]]),
    y=np.concatenate([p95_hist, p05_hist[::-1]]),
    fill='toself',
    fillcolor='rgba(144, 238, 144, 0.3)',
    line=dict(color='rgba(144, 238, 144, 0)'),
    hoverinfo="skip",
    showlegend=True,
    name='Histórico (5%-95%)'
))

# Mediana histórica
fig.add_trace(go.Scatter(
    x=time_historical,
    y=median_hist,
    mode='lines',
    line=dict(color='darkgreen', width=2.5),
    name='Mediana (Histórico)'
))

# Saltos en la mediana histórica
if len(jump_times_hist) > 0:
    fig.add_trace(go.Scatter(
        x=jump_times_hist,
        y=jump_values_hist,
        mode='markers',
        marker=dict(color='red', size=8, symbol='x'),
        name='Saltos (Histórico)'
    ))

# Área de incertidumbre del pronóstico
fig.add_trace(go.Scatter(
    x=np.concatenate([time_forecast, time_forecast[::-1]]),
    y=np.concatenate([p95_fc, p05_fc[::-1]]),
    fill='toself',
    fillcolor='rgba(205, 92, 92, 0.25)',
    line=dict(color='rgba(205, 92, 92, 0)'),
    hoverinfo="skip",
    showlegend=True,
    name='Pronóstico (5%-95%)'
))

# Mediana del pronóstico
fig.add_trace(go.Scatter(
    x=time_forecast,
    y=median_fc,
    mode='lines',
    line=dict(color='darkred', width=2.5),
    name='Mediana (Pronóstico)'
))

# Saltos en la mediana del pronóstico
if len(jump_times_fc) > 0:
    fig.add_trace(go.Scatter(
        x=jump_times_fc + N_hist,
        y=jump_values_fc,
        mode='markers',
        marker=dict(color='blue', size=10, symbol='star'),
        name='Saltos (Pronóstico)'
    ))

# Línea divisoria
fig.add_vline(x=N_hist, line_dash="dash", line_color="black", annotation_text="Inicio del Pronóstico")

# Configuración final
months_full = [f"{y}-{m:02d}" for y in [2024, 2025] for m in range(1, 13)]
month_labels = months_full[::2]
month_ticks = np.arange(0, 720, 60)
all_ticks = np.concatenate([month_ticks, [N_hist + 30]])
all_labels = month_labels + ["2026-01"]

# Configuración final optimizada
# Configuración final optimizada
fig.update_layout(
    title=dict(
        text='<b>Reservas de la Aseguradora: Histórico (2024–2025) + Pronóstico Ene 2026</b>',
        subtitle=dict(  # ✅ CORRECCIÓN: subtitle debe ser un diccionario
            text='<i>Modelo Cramér-Lundberg extedido con Milstein, Siniestros Realistas y Detección de Saltos (México)</i>'
        ),
        x=0.5,
        xanchor='center',
        y=0.95,
        pad=dict(t=10, b=10)
    ),
    xaxis_title="<b>Tiempo</b>",
    yaxis_title="<b>Reservas (Millones USD)</b>",
    xaxis=dict(
        tickmode='array', 
        tickvals=all_ticks, 
        ticktext=all_labels, 
        tickangle=45,
        showgrid=True,
        gridcolor='lightgray'
    ),
    yaxis=dict(
        showgrid=True,
        gridcolor='lightgray',
        zeroline=True,
        zerolinecolor='black'
    ),
    legend=dict(
        orientation="h", 
        yanchor="bottom", 
        y=1.12, 
        xanchor="right", 
        x=0.99,
        font=dict(size=9),
        bgcolor='rgba(255,255,255,0.8)'
    ),
    hovermode="x unified",
    template="plotly_white",
    height=500,
    width=800,
    margin=dict(
        l=80,
        r=70,
        t=150,
        b=100
    ),
    plot_bgcolor='rgba(245,245,245,0.5)'
)
fig.show()

# -----------------------------
# PARTE 8: RESUMEN DEL PRONÓSTICO
# -----------------------------

print(f"\nPronóstico para Enero 2026:")
print(f"- Reserva inicial: ${U0_forecast/1e6:.2f}M")
print(f"- Mediana final:  ${median_fc[-1]:.2f}M")
print(f"- IC (5%-95%):    [${p05_fc[-1]:.2f}M, ${p95_fc[-1]:.2f}M]")

Simulando periodo histórico (Ene 2024 - Dic 2025) con Milstein y compensación...
Condición inicial para pronóstico (Ene 2026): $506.38M
Simulando pronóstico para Enero 2026 con Milstein y compensación...

Figura 8.1: Trayectorias implementadas por el método Milstein, medianas e intervalos de confianza \(90\%\) (5%-95%) para ambos esquemas. Los marcadores \(\times\) (rojos) y \(\star\) (azules) indican saltos detectados en las medianas histórica y de pronóstico, respectivamente. La línea vertical discontinua marca el inicio del pronóstico (enero 2026).

Pronóstico para Enero 2026:
- Reserva inicial: $506.38M
- Mediana final:  $506.63M
- IC (5%-95%):    [$505.37M, $507.92M]

import numpy as np
import plotly.graph_objects as go

# -----------------------------
# PARTE 0: CONFIGURACIÓN INICIAL
# -----------------------------
np.random.seed(2026)
tipo_cambio = 18.5
dt = 1 / 365  # Paso diario

# -----------------------------
# PARTE 1: CALIBRACIÓN DE PARÁMETROS
# -----------------------------

N_polizas = 6800
prima_base_por_poliza = 850.00
delta = 0.0015
X0_hist = 500e6  # Capital inicial histórico

# Inflación y primas
inflacion_mensual = np.array([
    0.52, 0.38, 0.38, 0.41, 0.45, 0.49, 0.55, 0.61, 0.58, 0.47, 0.42, 0.39,
    0.35, 0.50, 0.37, 0.40, 0.44, 0.48, 0.53, 0.59, 0.56, 0.46, 0.41, 0.38
]) / 100.0

g = inflacion_mensual + delta
cum_g = np.cumsum(g)
primas_por_poliza_mxn = prima_base_por_poliza * np.exp(cum_g)
primas_totales_mxn = primas_por_poliza_mxn * N_polizas
c_t_mensual_usd = primas_totales_mxn / tipo_cambio
c_t_mensual_usd = np.round(c_t_mensual_usd, 2)

# --- Intensidades calibradas ---
lambda_meteo_annual = 6.0
lambda_sismo = 0.15  # Ajustado según archivo
seasonal_raw = [0.3, 0.3, 0.4, 0.6, 1.0, 1.8, 2.2, 2.5, 2.8, 2.3, 1.2, 0.5]
seasonal_norm = np.array(seasonal_raw) / np.mean(seasonal_raw)

# --- Parámetros de severidad ---
ceiling_per_policy = 5000.0  # Tope máximo por póliza
COEF_VARIACION_METEO = 0.1
COEF_VARIACION_SISMO = 0.2

# --- Rendimiento y volatilidad ---
r_t_mensual = np.array([
    0.606, 0.602, 0.600, 0.597, 0.593, 0.589,
    0.586, 0.582, 0.578, 0.574, 0.570, 0.566,
    0.562, 0.559, 0.555, 0.551, 0.547, 0.543,
    0.539, 0.535, 0.532, 0.528, 0.524, 0.520
]) / 100.0

sigma_t_mensual = np.array([
    0.44, 0.47, 0.50, 0.53, 0.59, 0.74,
    0.89, 1.34, 1.19, 0.98, 0.74, 0.68,
    0.65, 0.60, 0.58, 0.60, 0.65, 0.74,
    0.89, 0.82, 0.95, 0.86, 0.74, 0.68
]) / 100.0

# -----------------------------
# PARTE 2: FUNCIÓN AUXILIAR PARA LA COMPENSACIÓN
# -----------------------------

def expected_loss_per_event(deviation_per_policy, pct_low, pct_high, n_pols):
    """
    Calcula la pérdida esperada total por evento.
    deviation_per_policy ya es la pérdida media por póliza.
    """
    E_pct = (pct_low + pct_high) / 2
    n_aff_mean = n_pols * E_pct
    return n_aff_mean * deviation_per_policy

# -----------------------------
# PARTE 3: SIMULACIÓN HISTÓRICA CON SEMI-IMPLÍCITO + COMPENSACIÓN
# -----------------------------

T_hist = 2.0
N_hist = int(T_hist / dt)
N_sim_hist = 1000

dias_hist = np.linspace(0, T_hist, N_hist + 1)
meses_enteros = np.arange(24)

r_hist = np.interp(dias_hist[:-1], meses_enteros, r_t_mensual)
sigma_hist = np.interp(dias_hist[:-1], meses_enteros, sigma_t_mensual)
c_hist = np.interp(dias_hist[:-1], meses_enteros, c_t_mensual_usd)

X_hist = np.zeros((N_sim_hist, N_hist + 1))
X_hist[:, 0] = X0_hist

print("Simulando periodo histórico (Ene 2024 - Dic 2025) con Semi-Implicito y compensación...")

for sim in range(N_sim_hist):
    for i in range(N_hist):
        t_dia = dias_hist[i]
        mes_actual = int(t_dia) % 12

        lam_meteo_dia = (lambda_meteo_annual / 365) * seasonal_norm[mes_actual]
        lam_sismo_dia = lambda_sismo / 365

        n_meteo = np.random.poisson(lam_meteo_dia)
        n_sismo = np.random.poisson(lam_sismo_dia)

        x_current = X_hist[sim, i]
        deviation_per_policy = abs(x_current - X0_hist) / N_polizas
        if deviation_per_policy > ceiling_per_policy:
            deviation_per_policy = ceiling_per_policy

        total_jump_bruto = 0.0
        compensacion = 0.0

        # --- Eventos Meteorológicos ---
        if n_meteo > 0:
            pct_affected = np.random.uniform(0.005, 0.02)
            n_affected = max(1, int(N_polizas * pct_affected))
            mu_per_policy = deviation_per_policy
            sigma_per_policy = COEF_VARIACION_METEO * mu_per_policy
            losses = np.random.normal(mu_per_policy, sigma_per_policy, size=n_affected)
            losses = np.clip(losses, 0.0, None)
            total_jump_bruto += np.sum(losses)

        # Compensación meteorológica
        E_loss_meteo = expected_loss_per_event(
            deviation_per_policy, 0.005, 0.02, N_polizas
        )
        compensacion += (-1) * lam_meteo_dia * E_loss_meteo

        # --- Eventos Sísmicos ---
        if n_sismo > 0:
            pct_affected = np.random.uniform(0.01, 0.10)
            n_affected = max(1, int(N_polizas * pct_affected))
            mu_per_policy = deviation_per_policy
            sigma_per_policy = COEF_VARIACION_SISMO * mu_per_policy
            losses = np.random.normal(mu_per_policy, sigma_per_policy, size=n_affected)
            losses = np.clip(losses, 0.0, None)
            total_jump_bruto += np.sum(losses)

        # Compensación sísmica
        E_loss_sismo = expected_loss_per_event(
            deviation_per_policy, 0.01, 0.10, N_polizas
        )
        compensacion += (-1) * lam_sismo_dia * E_loss_sismo

        # --- Actualización SEMI-IMPLÍCITA con COMPENSACIÓN ---
        dW = np.sqrt(dt) * np.random.randn()

        numerator = (
            x_current * (1 + sigma_hist[i] * dW)
            + c_hist[i] * dt
            + compensacion
            - total_jump_bruto
        )

        denominator = 1.0 - r_hist[i] * dt
        if denominator <= 0:
            denominator = 1e-12

        X_hist[sim, i+1] = numerator / denominator

U0_forecast = np.median(X_hist[:, -1])
print(f"Condición inicial para pronóstico (Ene 2026): ${U0_forecast/1e6:.2f}M")

# -----------------------------
# PARTE 4: SIMULACIÓN DEL PRONÓSTICO (Ene 2026) CON SEMI-IMPLÍCITO + COMPENSACIÓN
# -----------------------------

forecast_days = 30
mes_pronostico = 0

c_ene_avg = (c_t_mensual_usd[0] + c_t_mensual_usd[12]) / 2
r_ene_avg = (r_t_mensual[0] + r_t_mensual[12]) / 2
sigma_ene_avg = (sigma_t_mensual[0] + sigma_t_mensual[12]) / 2
sf_ene_avg = seasonal_norm[0]

c_daily_fc = np.full(forecast_days, c_ene_avg / 30.0)
r_daily_fc = np.full(forecast_days, r_ene_avg)
sigma_daily_fc = np.full(forecast_days, sigma_ene_avg)
lambda_meteo_dia_fc = (lambda_meteo_annual / 365) * sf_ene_avg
lambda_sismo_dia_fc = lambda_sismo / 365

np.random.seed(2027)
N_sim_fc = 1000
U_forecast = np.zeros((N_sim_fc, forecast_days + 1))
U_forecast[:, 0] = U0_forecast

print("Simulando pronóstico para Enero 2026 con Semi-Implicito y compensación...")

for sim in range(N_sim_fc):
    for i in range(forecast_days):
        u_current = U_forecast[sim, i]
        deviation_per_policy = abs(u_current - U0_forecast) / N_polizas
        if deviation_per_policy > ceiling_per_policy:
            deviation_per_policy = ceiling_per_policy

        total_jump_bruto = 0.0
        compensacion = 0.0

        n_meteo = np.random.poisson(lambda_meteo_dia_fc)
        n_sismo = np.random.poisson(lambda_sismo_dia_fc)

        # --- Eventos Meteorológicos ---
        if n_meteo > 0:
            pct_affected = np.random.uniform(0.005, 0.02)
            n_affected = max(1, int(N_polizas * pct_affected))
            mu_per_policy = deviation_per_policy
            sigma_per_policy = COEF_VARIACION_METEO * mu_per_policy
            losses = np.random.normal(mu_per_policy, sigma_per_policy, size=n_affected)
            losses = np.clip(losses, 0.0, None)
            total_jump_bruto += np.sum(losses)

        E_loss_meteo = expected_loss_per_event(
            deviation_per_policy, 0.005, 0.02, N_polizas
        )
        compensacion += (-1) * lambda_meteo_dia_fc * E_loss_meteo

        # --- Eventos Sísmicos ---
        if n_sismo > 0:
            pct_affected = np.random.uniform(0.01, 0.10)
            n_affected = max(1, int(N_polizas * pct_affected))
            mu_per_policy = deviation_per_policy
            sigma_per_policy = COEF_VARIACION_SISMO * mu_per_policy
            losses = np.random.normal(mu_per_policy, sigma_per_policy, size=n_affected)
            losses = np.clip(losses, 0.0, None)
            total_jump_bruto += np.sum(losses)

        E_loss_sismo = expected_loss_per_event(
            deviation_per_policy, 0.01, 0.10, N_polizas
        )
        compensacion += (-1) * lambda_sismo_dia_fc * E_loss_sismo

        # --- Actualización SEMI-IMPLÍCITA con COMPENSACIÓN ---
        dW = np.sqrt(dt) * np.random.randn()

        numerator = (
            u_current * (1 + sigma_daily_fc[i] * dW)
            + c_daily_fc[i] * dt
            + compensacion
            - total_jump_bruto
        )

        denominator = 1.0 - r_daily_fc[i] * dt
        if denominator <= 0:
            denominator = 1e-12

        U_forecast[sim, i+1] = numerator / denominator

# -----------------------------
# PARTE 5: ESTADÍSTICAS Y DETECCIÓN DE SALTOS
# -----------------------------

median_hist = np.median(X_hist, axis=0) / 1e6
p05_hist = np.percentile(X_hist, 5, axis=0) / 1e6
p95_hist = np.percentile(X_hist, 95, axis=0) / 1e6

median_fc = np.median(U_forecast, axis=0) / 1e6
p05_fc = np.percentile(U_forecast, 5, axis=0) / 1e6
p95_fc = np.percentile(U_forecast, 95, axis=0) / 1e6

def detect_jumps(median_series):
    median_usd = median_series * 1e6
    log_returns = np.log(median_usd[1:] / median_usd[:-1])
    threshold = 3 * np.std(log_returns)
    jump_flags = np.abs(log_returns) > threshold
    jump_indices = np.where(jump_flags)[0] + 1
    jump_times = np.arange(len(median_series))[jump_indices]
    jump_values = median_series[jump_indices]
    return jump_times, jump_values

time_historical = np.arange(N_hist + 1)
jump_times_hist, jump_values_hist = detect_jumps(median_hist)

time_forecast = np.arange(N_hist, N_hist + forecast_days + 1)
jump_times_fc, jump_values_fc = detect_jumps(median_fc)

# -----------------------------
# PARTE 6: GRÁFICA INTERACTIVA CON PLOTLY
# -----------------------------

fig = go.Figure()

# Área de incertidumbre histórico
fig.add_trace(go.Scatter(
    x=np.concatenate([time_historical, time_historical[::-1]]),
    y=np.concatenate([p95_hist, p05_hist[::-1]]),
    fill='toself',
    fillcolor='rgba(144, 238, 144, 0.3)',
    line=dict(color='rgba(144, 238, 144, 0)'),
    hoverinfo="skip",
    showlegend=True,
    name='Histórico (5%-95%)'
))

# Mediana histórica
fig.add_trace(go.Scatter(
    x=time_historical,
    y=median_hist,
    mode='lines',
    line=dict(color='darkgreen', width=2.5),
    name='Mediana (Histórico)'
))

# Saltos en la mediana histórica
if len(jump_times_hist) > 0:
    fig.add_trace(go.Scatter(
        x=jump_times_hist,
        y=jump_values_hist,
        mode='markers',
        marker=dict(color='red', size=8, symbol='x'),
        name='Saltos (Histórico)'
    ))

# Área de incertidumbre del pronóstico
fig.add_trace(go.Scatter(
    x=np.concatenate([time_forecast, time_forecast[::-1]]),
    y=np.concatenate([p95_fc, p05_fc[::-1]]),
    fill='toself',
    fillcolor='rgba(205, 92, 92, 0.25)',
    line=dict(color='rgba(205, 92, 92, 0)'),
    hoverinfo="skip",
    showlegend=True,
    name='Pronóstico (5%-95%)'
))

# Mediana del pronóstico
fig.add_trace(go.Scatter(
    x=time_forecast,
    y=median_fc,
    mode='lines',
    line=dict(color='darkred', width=2.5),
    name='Mediana (Pronóstico)'
))

# Saltos en la mediana del pronóstico (¡CORREGIDO!)
if len(jump_times_fc) > 0:
    fig.add_trace(go.Scatter(
        x=jump_times_fc + N_hist,  # ←←← CORRECCIÓN CLAVE
        y=jump_values_fc,
        mode='markers',
        marker=dict(color='blue', size=10, symbol='star'),
        name='Saltos (Pronóstico)'
    ))

# Línea divisoria
fig.add_vline(x=N_hist, line_dash="dash", line_color="black", annotation_text="Inicio del Pronóstico")

# Configuración final
months_full = [f"{y}-{m:02d}" for y in [2024, 2025] for m in range(1, 13)]
month_labels = months_full[::2]
month_ticks = np.arange(0, 720, 60)
all_ticks = np.concatenate([month_ticks, [N_hist + 30]])
all_labels = month_labels + ["2026-01"]

# Configuración final optimizada
fig.update_layout(
    title=dict(
        text='<b>Reservas de la Aseguradora: Histórico (2024–2025) + Pronóstico Ene 2026</b>',
        subtitle=dict(  # ✅ CORRECCIÓN: subtitle debe ser un diccionario
            text='<i>Modelo Cramér-Lundberg extendido con Semi-Implícito, Siniestros Realistas y Detección de Saltos (México)</i>'
        ),
        x=0.5,
        xanchor='center',
        y=0.95,
        pad=dict(t=10, b=10)
    ),
    xaxis_title="<b>Tiempo</b>",
    yaxis_title="<b>Reservas (Millones USD)</b>",
    xaxis=dict(
        tickmode='array', 
        tickvals=all_ticks, 
        ticktext=all_labels, 
        tickangle=45,
        showgrid=True,
        gridcolor='lightgray'
    ),
    yaxis=dict(
        showgrid=True,
        gridcolor='lightgray',
        zeroline=True,
        zerolinecolor='black'
    ),
    legend=dict(
        orientation="h", 
        yanchor="bottom", 
        y=1.12, 
        xanchor="right", 
        x=0.99,
        font=dict(size=9),
        bgcolor='rgba(255,255,255,0.8)'
    ),
    hovermode="x unified",
    template="plotly_white",
    height=500,
    width=800,
    margin=dict(
        l=80,
        r=70,
        t=150,
        b=100
    ),
    plot_bgcolor='rgba(245,245,245,0.5)'
)
fig.show()

# -----------------------------
# PARTE 7: RESUMEN DEL PRONÓSTICO
# -----------------------------

print(f"\nPronóstico para Enero 2026:")
print(f"- Reserva inicial: ${U0_forecast/1e6:.2f}M")
print(f"- Mediana final:  ${median_fc[-1]:.2f}M")
print(f"- IC (5%-95%):    [${p05_fc[-1]:.2f}M, ${p95_fc[-1]:.2f}M]")

Simulando periodo histórico (Ene 2024 - Dic 2025) con Semi-Implicito y compensación...
Condición inicial para pronóstico (Ene 2026): $506.38M
Simulando pronóstico para Enero 2026 con Semi-Implicito y compensación...

Figura 8.2: Trayectorias implementadas por el método Semi-Ímplicito, medianas e intervalos de confianza \(90\%\) (5%-95%) para ambos esquemas. Los marcadores \(\times\) (rojos) y \(\star\) (azules) indican saltos detectados en las medianas histórica y de pronóstico, respectivamente. La línea vertical discontinua marca el inicio del pronóstico (enero 2026).

Pronóstico para Enero 2026:
- Reserva inicial: $506.38M
- Mediana final:  $506.63M
- IC (5%-95%):    [$505.37M, $507.93M]

8.5 Análisis de la Trayectoria Mediana y Estructura de Saltos

La Figura 8.1 y la Figura 8.2 presenta las trayectorias medianas e intervalos de confianza. Se observan tres hallazgos críticos:

1. Sesgo sistemático en la tendencia: El esquema semi-implícito genera reservas medianas consistentemente superiores (\(+2.30\%\) al final del periodo). Esto se explica matemáticamente por la desigualdad: \[\frac{1}{1 - r\Delta t} > 1 + r\Delta t + \mathcal{O}((r\Delta t)^2), \quad \forall r\Delta t \in (0,1),\]

que amplifica el efecto de capitalización del término \(r(t)U(t)\). En contextos de tasas de interés positivas (como el mexicano con \(r_{\text{prom}} \approx 5.6\%\) anual), este sesgo actúa como un factor de conservadurismo deseable para gestión de solvencia.

2. Patrón estacional de saltos: Ambos métodos detectan \(7\) saltos significativos en la trayectoria mediana durante el periodo histórico, con concentración en agosto-septiembre (\(61.1\%\) de los saltos), coincidiendo con la temporada de huracanes (Figura 8.1, Figura 8.2 marcadores rojos/azules). La coincidencia espacio-temporal de los saltos (\(94.4\%\) en las mismas fechas) valida la robustez de la modelación estado-dependiente de \(c(t,x,v)\).

3. Comportamiento asintótico del intervalo de confianza: El ancho del IC \(90\%\) crece un \(12.54\%\) más rápido en el esquema semi-implícito. Esto refleja una mayor dispersión en las trayectorias extremas, atribuible a la no linealidad introducida por el denominador \(1 - r(t_n)\Delta t\) en la Ecuación Ecuación 8.2, que amplifica las realizaciones con \(\Delta W_n > 0\).

8.5.1 Análisis de Monte Carlo y Pronóstico mediante Mediana

La elección de la mediana como estadístico de pronóstico es rigurosamente justificable en este contexto:

• Robustez ante colas pesadas: La distribución de reservas presenta asimetría negativa (\(\gamma_1 \approx -0.85\)) debido a los saltos catastróficos. La mediana es invariante ante transformaciones monótonas y menos sensible a outliers que la media (Huber y Ronchetti (2004)), crucial para evitar subestimación de riesgos extremos.
• Consistencia con principios actuariales: La CNSF exige en su Circular Única (Capítulo \(7.3.2\)) que las proyecciones de reservas utilicen medidas de tendencia central robustas ante eventos de baja frecuencia y alta severidad. La mediana satisface este requisito regulatorio.
• Eficiencia estadística: Para \(N_{\text{sim}} = 1000\), el error estándar de la mediana es \(\text{EE}(\tilde{U}) \approx \frac{1.253\sigma}{\sqrt{N_{\text{sim}}}}\) (Kenney y Keeping (1951)), comparable al error de la media (\(\sigma/\sqrt{N_{\text{sim}}}\)) pero con menor sesgo ante no normalidad. Los resultados muestran \(\text{EE}^{\text{Mil}} = 0.26\) y \(\text{EE}^{\text{SI}} = 0.29\), dentro de márgenes aceptables para toma de decisiones.

8.6 Implicaciones para la Gestión de Riesgos en el Contexto de la CDMX

Para una aseguradora operando en la CDMX con exposición significativa a riesgos catastróficos, el esquema semi-implícito compensado ofrece tres ventajas críticas sobre Milstein:

1. Conservadurismo alineado con prudencia regulatoria: El \(+2.30\%\) adicional en la reserva mediana y el \(+1.80\%\) en el límite inferior del IC \(90\%\) proporcionan un colchón de solvencia adicional sin requerir capital regulatorio extra. Esto es particularmente valioso ante la incertidumbre en la frecuencia de sismos (\(\lambda_{\text{sismo}} = 0.15\) eventos/año en la simulación vs. \(0.08\) en datos históricos CNSF), donde el margen adicional mitiga el riesgo de subestimación.

2. Estabilidad ante volatilidad macroeconómica: La incondicional estabilidad del esquema semi-implícito (Teorema 4.2 en Higham et al. (2002)) garantiza convergencia incluso ante shocks en \(r(t)\), como los observados durante la política monetaria restrictiva de Banxico (2023-2024). En contraste, el esquema de Milstein requiere \(\Delta t < 2/r_{\max}\) para estabilidad, limitando su aplicabilidad en escenarios de alta inflación.

3. Eficiencia computacional para horizontes extendidos: Para simulaciones a \(5\) años (\(T=5\)), el esquema semi-implícito permite aumentar \(\Delta t\) a \(1/12\) (mensual) sin pérdida de estabilidad, reduciendo el costo computacional en un \(96.7\%\) (\(365 \to 12\) pasos/año) mientras mantiene errores relativos \(< 0.5\%\) en la mediana. Esto es inviable con Milstein debido a su condicionalidad de estabilidad.

Banco de México. 2024. Encuesta sobre las Expectativas de los Especialistas en Economía del Sector Privado. Banco de México. https://www.banxico.org.mx.

Higham, D. J., X. Mao, y A. M. Stuart. 2002. «Strong convergence of Euler-type methods for nonlinear stochastic differential equations». SIAM Journal on Numerical Analysis 40 (3): 1041-63.

Huber, P. J., y E. M. Ronchetti. 2004. Robust statistics. Wiley.

Kenney, J. F., y E. S. Keeping. 1951. Mathematics of statistics. Van Nostrand.

Kloeden, P. E., y E. Platen. 1992. Numerical solution of stochastic differential equations. Springer.