# -*- coding: utf-8 -*-
import os
import json
import warnings
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
from scipy.optimize import minimize
from scipy.stats import truncnorm, gaussian_kde
warnings.filterwarnings("ignore")
# =============================================================================
# CONFIGURACIÓN REPRODUCIBLE
# =============================================================================
SEED = 2026
rng = np.random.default_rng(SEED)
N_SIM = 1000
T_HIST = 2.0 # años
DT = 1.0 / 365.0 # paso diario expresado en años
N_HIST = int(round(T_HIST / DT))
DIAS_POR_ANIO = 365
TIPO_CAMBIO = 18.5
B_BOOT = 1000 # bootstrap por trayectoria
ALPHA_VALUES = np.arange(2.0, 5.01, 0.25)
VARPI_VALUES = [0.40, 0.45, 0.49]
ALPHA_REFERENCIA = 3.5
VARPI_REFERENCIA = 0.45
MIN_DETECCIONES = 30
OUTPUT_DIR = Path.cwd() / "salidas_milstein_agregado"
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)
# =============================================================================
# 1. CALIBRACIÓN DEL MODELO GENERADOR
# =============================================================================
N_POLIZAS = 6800
PRIMA_BASE_POR_POLIZA = 850.00
DELTA = 0.0015
X0_HIST = 500e6
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
primas_por_poliza_mxn = PRIMA_BASE_POR_POLIZA * np.exp(np.cumsum(g))
c_t_mensual_usd = np.round(primas_por_poliza_mxn * N_POLIZAS / TIPO_CAMBIO, 2)
# Dos mecanismos generadores del MODELO; no son clases inferidas.
lambda_1_annual = 6.0
lambda_2_annual = 0.15
seasonal_raw = np.array([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 = seasonal_raw / seasonal_raw.mean()
# Marcas generadoras del modelo. Sus etiquetas no pasan al bloque inferencial.
mu_1, sigma_1 = 12.5, 1.0
mu_2, sigma_2 = 1_742_000.0, 875_000.0
a_trunc_2 = (0.0 - mu_2) / sigma_2
EY_1 = np.exp(mu_1 + 0.5 * sigma_1**2)
EY_2 = truncnorm.mean(a_trunc_2, np.inf, loc=mu_2, scale=sigma_2)
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
t_years = np.arange(N_HIST + 1) * DT
t_months = t_years[:-1] * 12.0
months = np.arange(24)
r_hist = np.interp(t_months, months, r_t_mensual)
sigma_hist = np.interp(t_months, months, sigma_t_mensual)
c_hist = np.interp(t_months, months, c_t_mensual_usd)
# =============================================================================
# 2. SIMULACIÓN MONTE CARLO: 1000 TRAYECTORIAS
# =============================================================================
def simular_trayectorias():
"""
Genera X_hist. No devuelve un registro de tipos de evento.
Una vez generada la matriz, la inferencia usa únicamente X_hist.
"""
X = np.zeros((N_SIM, N_HIST + 1), dtype=float)
X[:, 0] = X0_HIST
for sim in range(N_SIM):
for i in range(N_HIST):
mes = int(t_months[i]) % 12
lam1_step = (lambda_1_annual / DIAS_POR_ANIO) * seasonal_norm[mes]
lam2_step = lambda_2_annual / DIAS_POR_ANIO
n1 = rng.poisson(lam1_step)
n2 = rng.poisson(lam2_step)
total_jump = 0.0
if n1 > 0:
total_jump += np.exp(rng.normal(mu_1, sigma_1, size=n1)).sum()
if n2 > 0:
total_jump += truncnorm.rvs(
a_trunc_2, np.inf, loc=mu_2, scale=sigma_2,
size=n2, random_state=rng
).sum()
# Compensación del salto por paso de un día.
compensacion_step = -(lam1_step * EY_1 + lam2_step * EY_2)
x = X[sim, i]
dW = np.sqrt(DT) * rng.normal()
drift = c_hist[i] * DT + r_hist[i] * x * DT
diffusion = sigma_hist[i] * x * dW
milstein_corr = 0.5 * sigma_hist[i]**2 * x * (dW**2 - DT)
X[sim, i + 1] = (
x + drift + compensacion_step + diffusion
+ milstein_corr - total_jump
)
return X
#print("=" * 84)
#print("SIMULACIÓN: Milstein Compensado")
#print("=" * 84)
X_hist = simular_trayectorias()
#print(f"Trayectorias simuladas: {X_hist.shape[0]}")
#print(f"Pasos por trayectoria: {X_hist.shape[1]-1}")
#print(f"Semilla: {SEED}")
# Desde aquí NO se usa ninguna etiqueta causal.
# La única entrada estadística es X_hist.
# =============================================================================
# 3. ESTIMACIÓN LOCAL DE VOLATILIDAD: BIPOWER VARIATION
# =============================================================================
def calcular_bipower_variation(trayectorias, K_window, dt):
M, N = trayectorias.shape
prev = trayectorias[:, :-1]
nxt = trayectorias[:, 1:]
if np.any(prev == 0):
raise ValueError("Se encontraron niveles U(t_{i-1})=0; no puede formarse el incremento relativo.")
r_tilde = np.zeros((M, N))
r_tilde[:, 1:] = (nxt - prev) / prev
abs_r = np.abs(r_tilde)
P = abs_r[:, :-1] * abs_r[:, 1:]
cumP = np.hstack([np.zeros((M, 1)), np.cumsum(P, axis=1)])
sigma_sq = np.zeros(N)
for n in range(K_window, N):
ventana = cumP[:, n] - cumP[:, n - K_window]
sigma_sq[n] = (np.pi / 2.0) * ventana.sum() / (K_window * M * dt)
sigma_sq[:K_window] = sigma_sq[K_window]
sigma_sq = np.maximum(sigma_sq, 1e-18)
return np.sqrt(sigma_sq), r_tilde
K_window = int(np.sqrt(N_HIST))
sigma_hat, r_tilde = calcular_bipower_variation(X_hist, K_window, DT)
# =============================================================================
# 4. DETECCIÓN POR UMBRAL Y EXTRACCIÓN DE SEVERIDADES
# =============================================================================
def construir_umbral(sigma_local, dt, alpha, varpi):
return alpha * sigma_local * (dt ** varpi)
def detectar_saltos_umbral(trayectorias, r_tilde, sigma_local, dt,
alpha, varpi, c_serie, r_serie):
M, N = trayectorias.shape
tau = construir_umbral(sigma_local, dt, alpha, varpi)
U_prev = trayectorias[:, :-1]
drift_rel = c_serie[np.newaxis, :N-1] / U_prev + r_serie[np.newaxis, :N-1]
e_tilde = np.zeros((M, N))
e_tilde[:, 1:] = r_tilde[:, 1:] - drift_rel * dt
D = np.zeros((M, N), dtype=bool)
D[:, 1:] = np.abs(e_tilde[:, 1:]) > tau[np.newaxis, 1:]
return D, e_tilde, tau
def extraer_severidades(trayectorias, e_tilde, D):
idx_traj, idx_t = np.where(D)
if len(idx_t) == 0:
return np.array([]), np.empty((0, 2), dtype=int), 0
U_prev = trayectorias[idx_traj, idx_t - 1]
Y_hat = -e_tilde[idx_traj, idx_t] * U_prev
# En el modelo los siniestros disminuyen la reserva: Y_hat debe ser > 0.
valid = np.isfinite(Y_hat) & (Y_hat > 0)
descartadas_signo = int((~valid).sum())
indices = np.column_stack([idx_traj[valid], idx_t[valid]])
return Y_hat[valid], indices, descartadas_signo
# =============================================================================
# 5. CALIBRACIÓN DEL DETECTOR SIN INFORMACIÓN
# =============================================================================
def tabla_calibracion(trayectorias, r_tilde, sigma_local):
rows = []
for varpi in VARPI_VALUES:
for alpha in ALPHA_VALUES:
D, e, tau = detectar_saltos_umbral(
trayectorias, r_tilde, sigma_local, DT,
float(alpha), float(varpi), c_hist, r_hist
)
sev, idx, fp_sign = extraer_severidades(trayectorias, e, D)
counts = np.bincount(idx[:, 0], minlength=N_SIM) if len(idx) else np.zeros(N_SIM, dtype=int)
z = tau[1:] / (sigma_local[1:] * np.sqrt(DT))
fp_gauss = float(np.mean(2.0 * stats.norm.sf(z)))
rows.append({
"alpha": float(alpha),
"varpi": float(varpi),
"n_detecciones_validas": int(len(sev)),
"detecciones_por_trayectoria": float(len(sev) / N_SIM),
"trayectorias_con_deteccion": int(np.sum(counts > 0)),
"mediana_severidad": float(np.median(sev)) if len(sev) else np.nan,
"media_severidad": float(np.mean(sev)) if len(sev) else np.nan,
"q90_severidad": float(np.quantile(sev, 0.90)) if len(sev) else np.nan,
"descartadas_por_signo": int(fp_sign),
"prob_fp_gauss_aprox_por_incremento": fp_gauss
})
return pd.DataFrame(rows)
def agregar_diagnosticos_estabilidad(tabla):
partes = []
for varpi, sub in tabla.groupby("varpi"):
sub = sub.sort_values("alpha").copy()
sub["cambio_rel_N_siguiente"] = (
sub["n_detecciones_validas"].shift(-1) - sub["n_detecciones_validas"]
).abs() / sub["n_detecciones_validas"].clip(lower=1)
sub["cambio_rel_mediana_siguiente"] = (
sub["mediana_severidad"].shift(-1) - sub["mediana_severidad"]
).abs() / sub["mediana_severidad"].abs().clip(lower=1.0)
partes.append(sub)
return pd.concat(partes, ignore_index=True)
calib = tabla_calibracion(X_hist, r_tilde, sigma_hat)
calib = agregar_diagnosticos_estabilidad(calib)
calib.to_csv(OUTPUT_DIR / "01_calibracion_umbral.csv", index=False)
alpha_opt, varpi_opt = ALPHA_REFERENCIA, VARPI_REFERENCIA
#print("\nCALIBRACIÓN DEL UMBRAL ")
#print("Se generó la tabla completa de sensibilidad (alpha, varpi).")
#print("No se fuerza un 'óptimo' sin información externa u objetivo de pérdida.")
#print(f"Valor de referencia para continuar: alpha={alpha_opt:.2f}, varpi={varpi_opt:.2f}")
#print("La elección final debe justificarse con estabilidad y sensibilidad.")
D_opt, e_opt, tau_opt = detectar_saltos_umbral(
X_hist, r_tilde, sigma_hat, DT,
alpha_opt, varpi_opt, c_hist, r_hist
)
severidades, indices, descartadas_signo = extraer_severidades(X_hist, e_opt, D_opt)
if len(severidades) < MIN_DETECCIONES:
raise RuntimeError("Muy pocas detecciones válidas para estimar una distribución.")
#print(f"Detecciones válidas: {len(severidades)}")
#print(f"Descartadas por signo/no finitas: {descartadas_signo}")
# =============================================================================
# 6. DISTRIBUCIÓN AGREGADA DE SEVERIDADES DETECTADAS
# =============================================================================
def ecdf(x):
xx = np.sort(np.asarray(x))
ff = np.arange(1, len(xx)+1) / len(xx)
return xx, ff
def kde_log_severidad(y):
"""
KDE sobre Z=log(Y). La densidad en la escala original se obtiene por:
f_Y(y) = f_Z(log y) / y.
"""
y = np.asarray(y)
y = y[np.isfinite(y) & (y > 0)]
z = np.log(y)
kde = gaussian_kde(z, bw_method="scott")
qlo, qhi = np.quantile(y, [0.001, 0.999])
grid_y = np.geomspace(max(qlo, np.min(y)), max(qhi, qlo * 1.001), 500)
dens_y = kde(np.log(grid_y)) / grid_y
return grid_y, dens_y, kde
def ajustar_familias_severidad_agregada(y):
"""
Ajustes PARAMÉTRICOS AGREGADOS opcionales.
Ninguna familia representa un tipo de siniestro.
"""
y = np.asarray(y)
y = y[(y > 0) & np.isfinite(y)]
fits = []
# Lognormal
sh, loc, sc = stats.lognorm.fit(y, floc=0)
ll = stats.lognorm.logpdf(y, sh, loc=0, scale=sc).sum()
fits.append(("Lognormal", 2, ll, {"sigma_log": sh, "mu_log": np.log(sc)}))
# Gamma
a, loc, sc = stats.gamma.fit(y, floc=0)
ll = stats.gamma.logpdf(y, a, loc=0, scale=sc).sum()
fits.append(("Gamma", 2, ll, {"shape": a, "scale": sc}))
# Weibull
c, loc, sc = stats.weibull_min.fit(y, floc=0)
ll = stats.weibull_min.logpdf(y, c, loc=0, scale=sc).sum()
fits.append(("Weibull", 2, ll, {"shape": c, "scale": sc}))
out = []
for name, k, ll, pars in fits:
out.append({
"familia": name,
"logLik": float(ll),
"AIC": float(2*k - 2*ll),
"parametros": json.dumps({k: float(v) for k, v in pars.items()})
})
return pd.DataFrame(out).sort_values("AIC")
sev_x, sev_F = ecdf(severidades)
grid_y, dens_y, kde_y = kde_log_severidad(severidades)
pd.DataFrame({"severidad_usd": severidades}).to_csv(
OUTPUT_DIR / "02_severidades_detectadas_agregadas.csv", index=False
)
pd.DataFrame({"y": sev_x, "F_empirica": sev_F}).to_csv(
OUTPUT_DIR / "03_ecdf_severidad.csv", index=False
)
pd.DataFrame({"y": grid_y, "densidad_kde": dens_y}).to_csv(
OUTPUT_DIR / "04_kde_severidad.csv", index=False
)
fits_sev = ajustar_familias_severidad_agregada(severidades)
fits_sev.to_csv(OUTPUT_DIR / "05_familias_severidad_agregada.csv", index=False)
# =============================================================================
# 7. BOOTSTRAP POR TRAYECTORIA PARA SEVERIDAD
# Condicional al detector calibrado.
# =============================================================================
sev_por_tray = [severidades[indices[:, 0] == j] for j in range(N_SIM)]
def bootstrap_severidad_por_trayectoria(B=B_BOOT):
boot = []
for b in range(B):
sampled = rng.integers(0, N_SIM, size=N_SIM)
partes = [sev_por_tray[j] for j in sampled if len(sev_por_tray[j]) > 0]
if not partes:
continue
s = np.concatenate(partes)
boot.append([
len(s),
np.mean(s),
np.median(s),
np.quantile(s, 0.25),
np.quantile(s, 0.75),
np.quantile(s, 0.90),
np.quantile(s, 0.95)
])
return np.asarray(boot)
boot_sev = bootstrap_severidad_por_trayectoria()
cols_boot = ["n", "media", "mediana", "q25", "q75", "q90", "q95"]
pd.DataFrame(boot_sev, columns=cols_boot).to_csv(
OUTPUT_DIR / "06_bootstrap_severidad.csv", index=False
)
def ic_percentil(v):
return np.quantile(v, [0.025, 0.975])
resumen_sev = {
"n_detectadas": int(len(severidades)),
"media": float(np.mean(severidades)),
"mediana": float(np.median(severidades)),
"q25": float(np.quantile(severidades, .25)),
"q75": float(np.quantile(severidades, .75)),
"q90": float(np.quantile(severidades, .90)),
"q95": float(np.quantile(severidades, .95)),
"IC95_media": ic_percentil(boot_sev[:, 1]).tolist(),
"IC95_mediana": ic_percentil(boot_sev[:, 2]).tolist()
}
# =============================================================================
# 8. TIEMPOS DE DETECCIÓN POR TRAYECTORIA
# =============================================================================
tiempos_por_tray = []
for j in range(N_SIM):
idx_j = indices[indices[:, 0] == j, 1]
# índice de malla -> días desde el inicio
t_j = np.sort(idx_j.astype(float))
tiempos_por_tray.append(t_j)
gaps_completos = []
gaps_por_tray = []
censurados = []
cens_por_tray = []
for j, tt in enumerate(tiempos_por_tray):
g = np.diff(tt) if len(tt) >= 2 else np.array([], dtype=float)
g = g[g > 0]
gaps_por_tray.append(g)
gaps_completos.extend(g.tolist())
# Censura terminal: válida para la espera iniciada en la última detección.
if len(tt) >= 1:
c = N_HIST - tt[-1]
if c >= 0:
censurados.append(float(c))
cens_por_tray.append(np.array([float(c)]))
else:
cens_por_tray.append(np.array([], dtype=float))
else:
# No se crea una falsa espera "entre saltos" para trayectorias sin detecciones.
cens_por_tray.append(np.array([], dtype=float))
gaps_completos = np.asarray(gaps_completos, dtype=float)
censurados = np.asarray(censurados, dtype=float)
pd.DataFrame({"gap_dias": gaps_completos}).to_csv(
OUTPUT_DIR / "07_gaps_completos.csv", index=False
)
pd.DataFrame({"censura_terminal_dias": censurados}).to_csv(
OUTPUT_DIR / "08_censuras_terminales.csv", index=False
)
# =============================================================================
# 9. KAPLAN-MEIER DESCRIPTIVO PARA LA DISTRIBUCIÓN MARGINAL DE GAPS
# =============================================================================
def kaplan_meier(completos, cens):
"""
Estimador KM aplicado a la colección marginal de gaps.
Los gaps múltiples de una trayectoria pueden ser dependientes; por ello se usa
como estimador descriptivo y la incertidumbre se obtiene por bootstrap de trayectoria.
"""
t = np.concatenate([np.asarray(completos), np.asarray(cens)])
event = np.concatenate([
np.ones(len(completos), dtype=int),
np.zeros(len(cens), dtype=int)
])
order = np.argsort(t)
t, event = t[order], event[order]
unique_times = np.unique(t)
S = 1.0
rows = [(0.0, 1.0, len(t), 0, 0)]
for u in unique_times:
at_risk = int(np.sum(t >= u))
d = int(np.sum((t == u) & (event == 1)))
c = int(np.sum((t == u) & (event == 0)))
if at_risk > 0 and d > 0:
S *= (1.0 - d / at_risk)
rows.append((float(u), float(S), at_risk, d, c))
return pd.DataFrame(rows, columns=["tiempo_dias", "S_KM", "en_riesgo", "eventos", "censuras"])
km = kaplan_meier(gaps_completos, censurados)
km.to_csv(OUTPUT_DIR / "09_kaplan_meier_gaps.csv", index=False)
# =============================================================================
# 10. AJUSTES PARAMÉTRICOS AGREGADOS CON CENSURA
# Exponencial / Weibull / Gamma. NO hay mezcla por tipo de evento.
# =============================================================================
def nll_exp(theta, comp, cens):
log_rate = theta[0]
rate = np.exp(log_rate)
ll = len(comp) * np.log(rate) - rate * np.sum(comp) - rate * np.sum(cens)
return -ll
def nll_weibull(theta, comp, cens):
log_k, log_s = theta
k, s = np.exp(log_k), np.exp(log_s)
ll = np.sum(stats.weibull_min.logpdf(comp, k, loc=0, scale=s))
ll += np.sum(stats.weibull_min.logsf(cens, k, loc=0, scale=s))
return -ll
def nll_gamma(theta, comp, cens):
log_k, log_s = theta
k, s = np.exp(log_k), np.exp(log_s)
ll = np.sum(stats.gamma.logpdf(comp, k, loc=0, scale=s))
ll += np.sum(stats.gamma.logsf(cens, k, loc=0, scale=s))
return -ll
def ajustar_tiempos_censurados(comp, cens):
if len(comp) < 10:
return pd.DataFrame(), {}
media = max(np.mean(comp), 1.0)
r1 = minimize(nll_exp, [np.log(1.0/media)], args=(comp, cens), method="Nelder-Mead")
rate = float(np.exp(r1.x[0]))
ll1 = float(-r1.fun)
r2 = minimize(nll_weibull, [0.0, np.log(media)], args=(comp, cens), method="Nelder-Mead")
k_w, s_w = map(float, np.exp(r2.x))
ll2 = float(-r2.fun)
r3 = minimize(nll_gamma, [0.0, np.log(media)], args=(comp, cens), method="Nelder-Mead")
k_g, s_g = map(float, np.exp(r3.x))
ll3 = float(-r3.fun)
rows = [
{"familia": "Exponencial", "k_param": 1, "logLik": ll1,
"AIC": 2 - 2*ll1, "parametros": json.dumps({"rate": rate})},
{"familia": "Weibull", "k_param": 2, "logLik": ll2,
"AIC": 4 - 2*ll2, "parametros": json.dumps({"shape": k_w, "scale": s_w})},
{"familia": "Gamma", "k_param": 2, "logLik": ll3,
"AIC": 4 - 2*ll3, "parametros": json.dumps({"shape": k_g, "scale": s_g})},
]
pars = {
"Exponencial": {"rate": rate},
"Weibull": {"shape": k_w, "scale": s_w},
"Gamma": {"shape": k_g, "scale": s_g}
}
return pd.DataFrame(rows).sort_values("AIC"), pars
fits_t, pars_t = ajustar_tiempos_censurados(gaps_completos, censurados)
fits_t.to_csv(OUTPUT_DIR / "10_familias_tiempos_censurados.csv", index=False)
# =============================================================================
# 11. BOOTSTRAP POR TRAYECTORIA PARA GAPS
# =============================================================================
def bootstrap_gaps_por_trayectoria(B=B_BOOT):
out = []
for b in range(B):
sampled = rng.integers(0, N_SIM, size=N_SIM)
gs = [gaps_por_tray[j] for j in sampled if len(gaps_por_tray[j]) > 0]
cs = [cens_por_tray[j] for j in sampled if len(cens_por_tray[j]) > 0]
g = np.concatenate(gs) if gs else np.array([])
c = np.concatenate(cs) if cs else np.array([])
if len(g) < 10:
continue
cv2 = np.var(g, ddof=1) / np.mean(g)**2 if len(g) > 1 else np.nan
out.append([
len(g), len(c), np.mean(g), np.median(g), cv2
])
return np.asarray(out)
boot_gaps = bootstrap_gaps_por_trayectoria()
pd.DataFrame(
boot_gaps,
columns=["n_completos", "n_censurados", "media_gap", "mediana_gap", "CV2"]
).to_csv(OUTPUT_DIR / "11_bootstrap_gaps.csv", index=False)
# =============================================================================
# 12. INTENSIDAD DE DETECCIONES: KERNEL CIRCULAR ANUAL
# =============================================================================
def kernel_intensidad_circular(doy_datos, grid_doy, bandwidth, M_norm, period=365):
"""
Intensidad de detecciones por trayectoria y por día.
Núcleo gaussiano circular, incluyendo 1/b.
"""
doy_datos = np.asarray(doy_datos, dtype=float)
out = np.zeros_like(grid_doy, dtype=float)
if len(doy_datos) == 0:
return out
for k, x in enumerate(grid_doy):
d = np.abs(doy_datos - x)
dc = np.minimum(d, period - d)
out[k] = np.sum(stats.norm.pdf(dc / bandwidth)) / (M_norm * bandwidth)
return out
all_detect_days = indices[:, 1].astype(float)
doy = np.mod(all_detect_days, DIAS_POR_ANIO)
grid_doy = np.arange(DIAS_POR_ANIO, dtype=float)
bandwidths = [30, 60, 91, 182]
intensidades = {}
for bw in bandwidths:
intensidades[bw] = kernel_intensidad_circular(doy, grid_doy, bw, N_SIM)
bandwidth_ref = 30
lambda_hat = intensidades[bandwidth_ref]
df_int = pd.DataFrame({"dia_del_anio": grid_doy})
for bw in bandwidths:
df_int[f"lambda_hat_bw_{bw}"] = intensidades[bw]
df_int.to_csv(OUTPUT_DIR / "12_intensidad_detectada_kernel.csv", index=False)
# =============================================================================
# 13. DIAGNÓSTICO DE CAMBIO DE TIEMPO
# =============================================================================
def intensidad_periodica_a_timeline(lambda_doy, total_days):
idx = np.arange(total_days + 1) % len(lambda_doy)
return lambda_doy[idx]
lambda_timeline = intensidad_periodica_a_timeline(lambda_hat, N_HIST)
Lambda_timeline = np.zeros_like(lambda_timeline, dtype=float)
Lambda_timeline[1:] = np.cumsum(lambda_timeline[:-1]) # integral discreta en días
rescaled = []
for tt in tiempos_por_tray:
if len(tt) >= 2:
s = Lambda_timeline[tt.astype(int)]
ds = np.diff(s)
rescaled.extend(ds[ds > 0].tolist())
rescaled = np.asarray(rescaled)
if len(rescaled):
u_rescaled = 1.0 - np.exp(-rescaled)
D_rescaling = float(stats.kstest(u_rescaled, "uniform").statistic)
else:
u_rescaled = np.array([])
D_rescaling = np.nan
pd.DataFrame({"u_rescaled": u_rescaled}).to_csv(
OUTPUT_DIR / "13_residuos_time_rescaling.csv", index=False
)
# =============================================================================
# 14. FIGURAS
# =============================================================================
# 14.1 Severidad: histograma + KDE agregada
fig, ax = plt.subplots(figsize=(9, 5.5))
ax.hist(severidades, bins=45, density=True, alpha=0.55, label="Severidades detectadas")
ax.plot(grid_y, dens_y, linewidth=2.0, label="KDE agregada")
ax.set_xlabel("Severidad detectada (USD)")
ax.set_ylabel("Densidad")
ax.set_title("Milstein Compensado: distribución agregada de severidades detectadas")
ax.legend()
fig.tight_layout()
fig.savefig(OUTPUT_DIR / "14_severidad_kde.png", dpi=180)
plt.close(fig)
# 14.2 ECDF
fig, ax = plt.subplots(figsize=(9, 5.5))
ax.step(sev_x, sev_F, where="post")
ax.set_xlabel("Severidad detectada (USD)")
ax.set_ylabel("F empírica")
ax.set_title("Milstein Compensado: función de distribución empírica de severidades")
fig.tight_layout()
fig.savefig(OUTPUT_DIR / "15_severidad_ecdf.png", dpi=180)
plt.close(fig)
# 14.3 Gaps completos: descriptivo
fig, ax = plt.subplots(figsize=(9, 5.5))
ax.hist(gaps_completos, bins=40, density=True, alpha=0.60)
ax.set_xlabel("Días entre detecciones consecutivas")
ax.set_ylabel("Densidad")
ax.set_title("Milstein Compensado: gaps completos (descriptivo)")
fig.tight_layout()
fig.savefig(OUTPUT_DIR / "16_gaps_histograma.png", dpi=180)
plt.close(fig)
# 14.4 Kaplan-Meier + supervivencias ajustadas
fig, ax = plt.subplots(figsize=(9, 5.5))
ax.step(km["tiempo_dias"], km["S_KM"], where="post", linewidth=2, label="Kaplan-Meier")
if pars_t:
xmax = max(
np.max(gaps_completos) if len(gaps_completos) else 1,
np.max(censurados) if len(censurados) else 1
)
xx = np.linspace(0, xmax, 500)
p = pars_t["Exponencial"]
ax.plot(xx, np.exp(-p["rate"] * xx), label="Exponencial (censurada)")
p = pars_t["Weibull"]
ax.plot(xx, stats.weibull_min.sf(xx, p["shape"], scale=p["scale"]),
label="Weibull (censurada)")
p = pars_t["Gamma"]
ax.plot(xx, stats.gamma.sf(xx, p["shape"], scale=p["scale"]),
label="Gamma (censurada)")
ax.set_xlabel("Tiempo entre detecciones (días)")
ax.set_ylabel("Supervivencia")
ax.set_ylim(0, 1.02)
ax.set_title("Milstein Compensado: supervivencia de los tiempos entre detecciones")
ax.legend()
fig.tight_layout()
fig.savefig(OUTPUT_DIR / "17_gaps_supervivencia.png", dpi=180)
plt.close(fig)
# 14.5 Intensidad detectada
fig, ax = plt.subplots(figsize=(9, 5.5))
for bw in bandwidths:
ax.plot(grid_doy, intensidades[bw], label=f"b={bw} días")
ax.set_xlabel("Día del año")
ax.set_ylabel("Intensidad estimada de detecciones por trayectoria y día")
ax.set_title("Milstein Compensado: sensibilidad del suavizado de intensidad")
ax.legend()
fig.tight_layout()
fig.savefig(OUTPUT_DIR / "18_intensidad_kernel.png", dpi=180)
plt.close(fig)
cv2 = (
float(np.var(gaps_completos, ddof=1) / np.mean(gaps_completos)**2)
if len(gaps_completos) > 1 else np.nan
)
resumen = {
"metodo": "Milstein Compensado",
"seed": SEED,
"N_sim": N_SIM,
"T_years": T_HIST,
"dt_years": DT,
"K_bipower": K_window,
"alpha": alpha_opt,
"varpi": varpi_opt,
"n_detecciones_validas": int(len(severidades)),
"severidad_media": resumen_sev["media"],
"severidad_mediana": resumen_sev["mediana"],
"IC95_media": resumen_sev["IC95_media"],
"IC95_mediana": resumen_sev["IC95_mediana"],
"n_gaps_completos": int(len(gaps_completos)),
"n_censuras_terminales": int(len(censurados)),
"media_gap_dias": float(np.mean(gaps_completos)) if len(gaps_completos) else None,
"mediana_gap_dias": float(np.median(gaps_completos)) if len(gaps_completos) else None,
"CV2_gap_descriptivo": cv2,
"familia_severidad_AIC_min": (
str(fits_sev.iloc[0]["familia"]) if len(fits_sev) else None
),
"familia_tiempo_AIC_min": (
str(fits_t.iloc[0]["familia"]) if len(fits_t) else None
),
"D_time_rescaling_diagnostico": D_rescaling
}
with open(OUTPUT_DIR / "19_resumen.json", "w", encoding="utf-8") as f:
json.dump(resumen, f, indent=2, ensure_ascii=False)
#print("\n" + "=" * 84)
#print("RESUMEN FINAL")
#print("=" * 84)
#for k, v in resumen.items():
# print(f"{k}: {v}")
#print(f"\nResultados guardados en: {OUTPUT_DIR.resolve()}")
