Empezamos por definir algunas nociones básicas sobre procesos de Hawkes:
Definición Un proceso de Hawkes \(\{N(t)\}_{t\geq0}\) definido por \(N(t)=\sum_{i<t}T_i\) y con historia \((\mathcal{H}_t:t \geq 0)\) es un proceso de Poisson doblemente estocástico, determinado por la función de intensidad condicional \[ \lambda(t|\mathcal{H}_t)=\lambda_0(t)+\sum_{i:t>t_i}\phi(t-t_i). \] Donde cada \(t_i\) es una realización de la variable aleatoria \(T_i\).
En general, tomamos \(\lambda(t):=\lambda(t|\mathcal{H}_t)\) salvo que se indique lo contrario. Al ser un proceso de Poisson doblemente estocástico, si tenemos la historia del proceso hasta el tiempo \(t\) se cumple con el siguiente sistema de ecuaciones: \[\mathbb{P}(N(t+h)-N(t) = m \mid \mathcal H_t) = \begin{cases} 1 - \lambda(t)\,h + o(h), & m = 0,\\[6pt] \lambda(t)\,h + o(h), & m = 1,\\[6pt] o(h), & m \ge 2, \end{cases}\]
Ahora, sea \(\Lambda(t)=\int_{0}^t\lambda(s)ds\) (compensador de la intensidad), entonces tenemos el siguiente importante resultado:
Teorema La función de verosimilitud del proceso de Hawkes \(N(*)\) con intensidad \(\lambda(t)\) en el tiempo \([0,t]\) va dada por \[ L=f(t_1,t_2,\ldots,t_{n(T)})=\prod_{i=1}^{n(T)}f(t_i|\mathcal H_{t_i} )=\prod_{i=1}^{n(T)}\lambda(t_i)e^{-\Lambda(t)} \] Otro resultado que debemos tener en cuenta es el reconocido Teorema de Bayes:
Teorema (Bayes) \[p(\theta|y)=\frac{p(y| \theta)p(\theta)}{\int_{\Theta}p(y|\theta)p(\theta)d\theta} \propto p(y| \theta)p(\theta)\] Donde \(p\) es una función de probabilidad.
Finalmente, utilizaremos el siguiente algoritmo del tipo Markov Chain Monte Carilo con el que estimaremos parámetros de un proceso de Hawkes.
Dado el estado actual \(\theta^{(b)}\) del parámetro \(\theta\), el siguiente estado \(\theta^{(b+1)}\) se obtiene mediante el siguiente procedimiento:
1. Generación de la propuesta
Simular un valor propuesto \(\theta' \sim
J(\theta \mid \theta^{(b)})\), donde \(J\) es una distribución de propuestas
arbitraria.
2. Cálculo de la tasa de aceptación
Evaluar la tasa de aceptación \(r\)
como
\[r = \frac{p(\theta' \mid y)}{p(\theta^{(b)} \mid y)} \cdot \frac{J(\theta^{(b)} \mid \theta')}{J(\theta' \mid \theta^{(b)})}.\]
Este ajuste incorpora el cociente de las probabilidades de la distribución de propuestas para corregir posibles asimetrías en \(J\).
3. Actualización del estado
Definir \(\theta^{(b+1)}\) de acuerdo
con la probabilidad de aceptación:
\[\theta^{(b+1)} = \begin{cases} \theta', & \text{con probabilidad } \min(1,r), \\ \theta^{(b)}, & \text{con probabilidad } 1 - \min(1,r). \end{cases}\]
Este algoritmo genera una cadena de Markov sobre el espacio de estados \(\Theta\) al que pertenece \(\theta\). En el Teorema 7.2 de [2] se prueba que \(p(\theta|y)\) es una distribución estacionaria de la cadena generada, y que en caso de ser una cadena irreducible entonces converge a dicha distribución.
Seguimos el procedimiento propuesto en [1].
Consideramos un proceso de Hawkes con función de intensidad (condicional) dada por \[ \lambda(t)=a+(\lambda_0-a)e^{-\beta t}+ \sum_{i>0:t_i<t}Y_ie^{-\beta(t-t_i)} \ , \ t>0\]
lambda <- function(t, a, lambda_0, beta, Y, P) {
suma_saltos <- sum(Y * exp(-beta * (t - P)) * (t > P))
resultado <- a + (lambda_0 - a) * exp(-beta * t) + suma_saltos
return(max(0, resultado)) }En el cuaderno anterior trabajamos con una intensidad muy similar, con la única diferencia de que tomábamos todo \(Y_i= \gamma\) constante. En cambio, ahora tomamos \(Y_i \sim \exp(1/b)\) la magnitud de salto dado por cada nuevo evento \(T_i=t_i\).
Por otro lado, no es difícil calcular la integral de la función de intensidad (condicional) en el intervalo \([0,T]\): \[\int_0^T\lambda(t)dt=aT+\frac{1}{\beta}(\lambda_0-a)(1-e^{-\beta T})+\frac{1}{\beta} \sum_{i=1}^{n(T)}Y_i(1-e^{-\beta (T-t_i)})\]
integral_lambda <- function(T_final, a, lambda_0, beta, Y, P) {
term_lineal <- a * T_final
term_exponencial <- ((lambda_0 - a) / beta) * (1 - exp(-beta * T_final))
term_saltos <- 0
for(i in 1:length(Y)) {
if(P[i] <= T_final) {
term_saltos <- term_saltos + Y[i] * (1 - exp(-beta * (T_final - P[i]))) / beta
}
}
return(term_lineal + term_exponencial + term_saltos)
}Objetivo Dado tiempos de eventos \(\{ t_i\}_{i=1}^{n(T)}\) en el intervalo \([0,T]\), queremos estimar el valor de los parámetros \((a,b,\beta)\)
Distribuciones previas Este es nuestro punto de partida para aplicar el Teorema de Bayes. Estas distribuciones se obtienen típicamente a partir del contexto del problema, opinión de expertos u otros estudios previos:
\(a \sim \text{Uniforme}(0,\lambda_0)\)
\(\beta \sim \text{Uniforme}(0,2)\)
\(b|\beta \sim \text{Uniforme}(0,\beta)\)
log_prior <- function(a, b, beta, lambda_0) {
# a ~ Uniform(0, lambda_0)
if (a <= 0 || a >= lambda_0) return(-Inf)
# beta ~ Uniform(0, 2)
if (beta <= 0 || beta >= 2) return(-Inf)
# b | beta ~ Uniform(0, beta)
if (b <= 0 || b >= beta) return(-Inf)
return(0)
}Ahora, sea \(\mathbf Y=(Y_i)_{i \geq 1}\) y \(\mathbf T=(T_i)_{i \geq 1}\), entonces tenemos que la función de verosimilitud del proceso va dada por \(L(\mathbf T)=p(\mathbf T | a,b, \beta, \mathbf{Y})\) y así \[p(a,b, \beta, \mathbf{Y} |\mathbf T ) \propto L(\mathbf T)p(a,b,\beta, \mathbf Y),\] aplicando el Teorema de Bayes.
Como \(\mathbf Y\) depende únicamente de \(b\), \(b\) depende únicamente de \(\beta\) y \(a\) y \(\beta\) son independientes de los demás parámetros entonces podemos escribir \[p(a,b,\beta, \mathbf Y)=p(\mathbf Y | b) p(b| \beta) p(a) p(b)\].
Remplazando con la expresión que conocemos para la función de verosimilitud, y las funciones de densidades de las distribuciones previas dadas para los parámetros, nos queda que \[ \begin{aligned} p(a,b, \beta, \mathbf{Y} |\mathbf T ) &\propto \prod_{i=0}^{n(T)} \lambda(t_i) e^{-\int_0^T \lambda_t dt} \times \prod_{i=0}^{n(T)} \frac{1}{b} e^{-\frac{Y_i}{b}} \cdot \frac{1}{\beta} \cdot \frac{1}{\lambda_0} \cdot \frac{1}{2} \\ &\propto \prod_{i=0}^{n(T)} \lambda(t_i) e^{-\int_0^T \lambda_t dt} \times \frac{1}{b^{n(T)}} e^{-\frac{1}{b} \sum_{i=1}^{n(T)} Y_i} \cdot \frac{1}{\beta}, \end{aligned} \]
siempre que \(a \in (0, \lambda_0)\), \(\beta \in (0,2)\) y \(b \in (0,\beta)\).
Y utilizando esta forma proporcional a $ p(a,b,,|T)$, calculamos la log-posterior
\[ \begin{aligned} \log (p(a, b, \beta, Y | T)) &= -n(T) \ln b - \frac{1}{b} \sum_{i=0}^{n(T)} Y_i + \sum_{i=0}^{n(T)} \ln(\lambda(t_i)) \\ &\quad - \int_0^T \lambda_t dt - \ln \beta + \text{const.}, \end{aligned}\]
log_posterior <- function(a, lambda_0, beta, Y, P, b, T_final) {
# 1. Verificar priors
log_prior_val <- log_prior(a, b, beta, lambda_0)
if (is.infinite(log_prior_val)) return(-Inf)
n_T <- length(P)
# 2. Prior de Y ~ Exp(b)
term_Y <- -n_T * log(b) - (1/b) * sum(Y)
# 3. Suma de log λ(P_i) - con verificación robusta
sum_log_lambda <- 0
for(i in 1:n_T) {
lambda_val <- lambda(P[i], a, lambda_0, beta, Y, P)
if(lambda_val <= 1e-10) return(-Inf) # Evitar underflow
sum_log_lambda <- sum_log_lambda + log(lambda_val)
}
# 4. Integral de la intensidad
integral_val <- integral_lambda(T_final, a, lambda_0, beta, Y, P)
if(is.na(integral_val) || integral_val < 0) return(-Inf)
term_integral <- -integral_val
# Combinar todo
log_post <- log_prior_val + term_Y + sum_log_lambda + term_integral
# Verificar resultado finito
if(is.na(log_post) || is.infinite(log_post)) return(-Inf)
return(log_post)
}Siguiendo el algoritmo de Metropolis-Hastings, las distribuciones \(J\) que se proponen van dadas por:
Distribuciones Propuestas Primero damos unas elecciones genéricas, las cuales más que por conocimiento del comportamiento de los parámetros son eficientes a la hora de aplicar el método:
\(a'|a \sim \text{N}(a,\sigma_a^2)\)
\(b'|b \sim \text{N}(b,\sigma_b^2)\)
\(\beta'| \beta \sim \text{N}(\beta,\sigma_\beta^2)\)
Y finalmente, damos una propuesta que se acomoda a nuestro conocimiento sobre las variables aleatorias \(Y_i\):
Luego, calculamos \[ \frac{J(a|a')}{J(a'|a)}=\frac{J(b|b')}{J(b'|b)}=\frac{J(\beta|\beta')}{J(\beta'|\beta)}=1 \ \text{(recuerde que la normal es simétrica)}\] y
\[ \frac{J(\mathbf Y|b)}{J(\mathbf Y'|b)}= \frac{\prod_{i=1}^{n(T)} \frac{1}{b} e^{-\frac{Y_i}{b}}}{\prod_{i=1}^{n(T)} \frac{1}{b} e^{-\frac{Y_i'}{b}}}= \exp\left( -\frac{1}{b} \sum_{i=1}^{n(T)} (Y_i - Y_i') \right)\].
De modo que los ratios de aceptación quedan \[ r_a=\frac{p(a',\mathbf{Y},b,\beta|\mathbf{T})}{p(a,\mathbf{Y},b,\beta|\mathbf{T})}\times1 \ , \ r_b=\frac{p(a,\mathbf{Y},b',\beta|\mathbf{T})}{p(a,\mathbf{Y},b,\beta|\mathbf{T})}\times1 \ , \ r_{\beta}=\frac{p(a,\mathbf{Y},b,\beta'|\mathbf{T})}{p(a,\mathbf{Y},b,\beta|\mathbf{T})}\times1\] y
\[r_{\mathbf{Y}}=\frac{p(a',\mathbf{Y},b,\beta|\mathbf{T})}{p(a,\mathbf{Y},b,\beta|\mathbf{T})}\times \exp\left( -\frac{1}{b} \sum_{i=1}^{n(T)} (Y_i - Y_i') \right) .\]
En la práctica se calculan estos ratios en escala \(\log\) , esto para tener mayor estabilidad númerica en la computación. Ahora veamos al algoritmo de Metropolis-Hastings que utilizaremos para nuestra estimación de los parámetros en el procesode Hawkes:
Dado tiempos de salto \(t = \{ t_0, t_1, ...,t_{n(T)} \}\). Inicialice \(a\), \(b\), \(\beta\) y \(\{ Y_i\}\) a su preferencia (por ejemplo, si tiene una idea de cuales pueden ser los valores reales), recordando que los valores \(Y_i\) provienen de una distribución exponencial con parámetro \(b^{-1}\). También tome \(k=0\), un número máximo de iteraciones \(M\) y escoja un valor de burn-in, esto es, todas las iteraciones antes del valor escogido de burn-in serán eliminados.
Es así como, mientras se cumpla \(k\) < \(M\) ejecute:
Muestreo de cada \(Y_i\):
Muestreo de \(a\):
Muestreo de \(b\):
Muestreo de \(\beta\):
Al finalizar, calcule el promedio de todos los pasos luego del burn-in para dar las estimaciones finales, esto es:
\(\hat{a}=\text{promedio}(a_{\text{burn-in}:M})\)
\(\hat{b}=\text{promedio}(b_{\text{burn-in}:M})\)
\(\hat{\beta}=\text{promedio}(\beta_{\text{burn-in}:M})\).
metropolis_hastings_hawkes <- function(P, n_iter, burn_in,
a_init, b_init, beta_init,
sigma_a, sigma_b, sigma_beta, T_final, lambda_0) {
# Inicializar cadenas
n_T <- length(P)
a_chain <- numeric(n_iter)
b_chain <- numeric(n_iter)
beta_chain <- numeric(n_iter)
# Estado actual
a_chain[1] <- a_init
b_chain[1] <- b_init
beta_chain[1] <- beta_init
Y_current <- rexp(n_T, rate = 1/b_init)
# Contadores de aceptación
aceptaciones <- list(a = 0, b = 0, beta = 0, Y = 0)
cat("Ejecutando MCMC...\n")
for(k in 2:n_iter) {
if(k %% 1000 == 0) cat("Iteración:", k, "/", n_iter, "\n")
# 1. ACTUALIZAR Y_i (variables latentes)
for(i in 1:n_T) {
Y_i_new <- rexp(1, rate = 1/b_chain[k-1])
Y_new <- Y_current
Y_new[i] <- Y_i_new
log_post_ratio <- log_posterior(a_chain[k-1], lambda_0, beta_chain[k-1], Y_new, P, b_chain[k-1], T_final) -
log_posterior(a_chain[k-1], lambda_0, beta_chain[k-1], Y_current, P, b_chain[k-1], T_final)
log_prop_ratio <- -(Y_current[i] - Y_i_new) / b_chain[k-1]
log_A_Yi <- log_post_ratio + log_prop_ratio
# Aceptar/rechazar de manera robusta
if(is.na(log_A_Yi) || log_A_Yi < -100) {
# Rechazar si hay problemas numéricos
} else if(log_A_Yi >= 0) {
Y_current[i] <- Y_i_new
aceptaciones$Y <- aceptaciones$Y + 1
} else {
A_Yi <- exp(log_A_Yi)
if(!is.na(A_Yi) && runif(1) < A_Yi) {
Y_current[i] <- Y_i_new
aceptaciones$Y <- aceptaciones$Y + 1
}
}
}
# 2. ACTUALIZAR a
a_new <- rnorm(1, a_chain[k-1], sigma_a)
if(a_new > 0 && a_new < lambda_0) {
log_A_a <- log_posterior(a_new, lambda_0, beta_chain[k-1], Y_current, P, b_chain[k-1], T_final) -
log_posterior(a_chain[k-1], lambda_0, beta_chain[k-1], Y_current, P, b_chain[k-1], T_final)
if(!is.na(log_A_a) && log_A_a >= 0) {
a_chain[k] <- a_new
aceptaciones$a <- aceptaciones$a + 1
} else if(!is.na(log_A_a) && log_A_a > -100) {
A_a <- exp(log_A_a)
if(!is.na(A_a) && runif(1) < A_a) {
a_chain[k] <- a_new
aceptaciones$a <- aceptaciones$a + 1
} else {
a_chain[k] <- a_chain[k-1]
}
} else {
a_chain[k] <- a_chain[k-1]
}
} else {
a_chain[k] <- a_chain[k-1]
}
# 3. ACTUALIZAR b
b_new <- rnorm(1, b_chain[k-1], sigma_b)
if(b_new > 0 && b_new < beta_chain[k-1]) {
log_A_b <- log_posterior(a_chain[k], lambda_0, beta_chain[k-1], Y_current, P, b_new, T_final) -
log_posterior(a_chain[k], lambda_0, beta_chain[k-1], Y_current, P, b_chain[k-1], T_final)
if(!is.na(log_A_b) && log_A_b >= 0) {
b_chain[k] <- b_new
aceptaciones$b <- aceptaciones$b + 1
} else if(!is.na(log_A_b) && log_A_b > -100) {
A_b <- exp(log_A_b)
if(!is.na(A_b) && runif(1) < A_b) {
b_chain[k] <- b_new
aceptaciones$b <- aceptaciones$b + 1
} else {
b_chain[k] <- b_chain[k-1]
}
} else {
b_chain[k] <- b_chain[k-1]
}
} else {
b_chain[k] <- b_chain[k-1]
}
# 4. ACTUALIZAR beta
beta_new <- rnorm(1, beta_chain[k-1], sigma_beta)
if(beta_new > 0 && beta_new < 2) {
log_A_beta <- log_posterior(a_chain[k], lambda_0, beta_new, Y_current, P, b_chain[k], T_final) -
log_posterior(a_chain[k], lambda_0, beta_chain[k-1], Y_current, P, b_chain[k], T_final)
if(!is.na(log_A_beta) && log_A_beta >= 0) {
beta_chain[k] <- beta_new
aceptaciones$beta <- aceptaciones$beta + 1
} else if(!is.na(log_A_beta) && log_A_beta > -100) {
A_beta <- exp(log_A_beta)
if(!is.na(A_beta) && runif(1) < A_beta) {
beta_chain[k] <- beta_new
aceptaciones$beta <- aceptaciones$beta + 1
} else {
beta_chain[k] <- beta_chain[k-1]
}
} else {
beta_chain[k] <- beta_chain[k-1]
}
} else {
beta_chain[k] <- beta_chain[k-1]
}
}
# Calcular estimaciones post burn-in
idx_post_burn <- (burn_in + 1):n_iter
a_hat <- mean(a_chain[idx_post_burn])
b_hat <- mean(b_chain[idx_post_burn])
beta_hat <- mean(beta_chain[idx_post_burn])
# Tasas de aceptación
tasa_aceptacion <- list(
a = aceptaciones$a / (n_iter - 1),
b = aceptaciones$b / (n_iter - 1),
beta = aceptaciones$beta / (n_iter - 1)
)
# RESULTADOS
cat("RESULTADOS MCMC - PROCESO DE HAWKES\n")
cat("ESTIMACIONES:\n")
cat(" a (tasa base): ", round(a_hat, 4), "\n")
cat(" b (escala saltos): ", round(b_hat, 4), "\n")
cat(" beta (decaimiento): ", round(beta_hat, 4), "\n")
cat("\nTASAS DE ACEPTACIÓN:\n")
cat(" a: ", round(tasa_aceptacion$a, 3), "\n")
cat(" b: ", round(tasa_aceptacion$b, 3), "\n")
cat(" beta: ", round(tasa_aceptacion$beta, 3), "\n")
return(list(
estimates = c(a = a_hat, b = b_hat, beta = beta_hat),
acceptance_rates = tasa_aceptacion,
chains = list(a = a_chain, b = b_chain, beta = beta_chain)
))
}Para poder experimentar con las estimaciones dadas por el algoritmo de Metropolis-Hastings, necesitamos también poder simular un proceso de Hawkes con kernel exponencial definido a partir de la función de intensidad que dimos antes en este documento. El siguiente algoritmo es tomado directamente de [3], y cumple justamente con dicha función:
simular_hawkes_exponencial <- function(N, lambda_inicial, a, alpha, beta, b) {
# Definimos vectores para guardar los valores generados por la simulación
T <- numeric(N + 1)
lambda_antes <- numeric(N + 1)
lambda_despues <- numeric(N + 1)
Y <- numeric(N + 1) # Nuevo vector para los saltos Y_i
# Condiciones iniciales
T[0] <- 0
lambda_antes[1] <- lambda_inicial
lambda_despues[1] <- lambda_inicial
Y[1] <- rexp(1, rate = 1/b) # Primer salto exponencial
for (i in 2:(N + 1)) {
# a) Tiempo hasta evento base (Poisson con tasa a)
u0 <- runif(1)
s0 <- -(1 / a) * log(u0)
# b) Tiempo hasta evento por autointensificación
u1 <- runif(1)
d <- 1 + (beta * log(u1)) / (lambda_despues[i - 1] - a)
# c) Seleccionar el tiempo mínimo
if (d > 0) {
s1 <- -(1 / beta) * log(d)
tau_i <- min(s0, s1)
} else {
tau_i <- s0
}
# d) Actualizar tiempo del evento
T[i] <- T[i - 1] + tau_i
# e) Intensidad justo antes del evento
lambda_antes[i] <- ((lambda_despues[i - 1] - a) * exp(-beta * tau_i)) + a
# f) Intensidad justo después del evento (salto exponencial)
Y[i] <- rexp(1, rate = 1/b) # Y_i ~ Exp(b⁻¹)
lambda_despues[i] <- lambda_antes[i] + Y[i]
}
return(list(
tiempos = T,
lambda_antes = lambda_antes,
lambda_despues = lambda_despues,
Y = Y
))
}Podemos hacer los experimentos para comprobar la efectividad del algoritmo en dos direcciones. Una primera dirección sería la siguiente:
Dado un conjunto de tiempos de los que no tenemos ninguna idea de donde provienen, pero que sospechamos que pueden ser ajustados dentro de un proceso de Hawkes con kernel exponencial. Entonces podemos utilizar nuestro algoritmo para encontrar los mejores parámetros (estimados) que hacen que los tiempos cuadren dentro de un proceso de Hakwes. Siendo así, tome los siguientes tiempos:
P <- c(0.5531861, 2.0379750, 3.8610338, 9.1488426, 9.4210556,
9.4344200, 10.8499658, 10.8574215, 19.4724278, 23.0319395)De la experiencia o de estudios anteriores, se dan valores iniciales para los parámetros a estimar, los cuales cumplan con las condiciones que usamos para formular las distribuciones previas.
Parámetros como el número de iteraciones que corremos el algoritmo y el de burn-in son difíciles de escoger de modo que se pueda asegurar cuándo son lo suficientemente buenos, pero con la experimentación podemos dar valores razonables. Podemos, por ejemplo, evaluar la distancia entre las estimaciones dadas entre iteraciones.
Finalmente, el valor del tiempo de observación \(T\) también tiene que ser escogido con cuidado, puesto que un valor demasiado cercano aumenta el riesgo de sobreestimar intensidad, y uno demasiado lejano en cambio aplica un sesgo hacia intensidades bajas y más cómputo. En [4], se recomienda tomar \(T=t_{n(T)}+c \frac{1}{\beta}\) con \(c\) entre 3 y 6.
# Parámetros razonables
lambda_0 <- 1.5
a_init <- 0.8
b_init <- 0.5
beta_init <- 1.2
T_final <- max(P) + 4
n_iter <- 12000
burn_in <- 4000Los valores \(\sigma\) se van ajustando a medida que experimentamos con los resulatdos del algoritmo. En la práctica, usamos la taza de aceptación de las nuevas propuestas que emite el algoritmo de Metropolis-Hastings, para ajustar nuestros valores \(\sigma\). Esto es, si la taza de aceptación es muy alta quiere decir que casi siempre tomamos la nueva propuesta y no estamos convergiendo a un valor, es así como debemos tomar un nuevo valor de \(\sigma\) para dar una mayor exploración del espacio de estados.
# Sigmas ajustados
sigma_a <- 0.3
sigma_b <- 0.3
sigma_beta <- 0.7Ejecutamos el algoritmo:
salida <- capture.output(
res <- metropolis_hastings_hawkes(
P = P, n_iter = n_iter, burn_in = burn_in,
a_init = a_init, b_init = b_init, beta_init = beta_init,
sigma_a = sigma_a, sigma_b = sigma_b, sigma_beta = sigma_beta,
T_final = T_final, lambda_0 = lambda_0
)
)
salida_filtrada <- salida[!grepl("Iteración", salida)]
cat(salida_filtrada, sep = "\n")## Ejecutando MCMC...
## RESULTADOS MCMC - PROCESO DE HAWKES
## ESTIMACIONES:
## a (tasa base): 0.2706
## b (escala saltos): 0.5393
## beta (decaimiento): 1.4089
##
## TASAS DE ACEPTACIÓN:
## a: 0.4
## b: 0.463
## beta: 0.504Ahora, queremos simular un proceso de Hawkes con los parámetros que estimamos vía el algoritmo MCMC. Utilizamos el algoritmo para simular procesos de Hawkes que presentamos anteriormente. Por la naturaleza aleatoria del proceso, no sería correcto basarnos en una sola simulación para comparar con los tiempos iniciales en el vector P. En cambio, tomamos 1000 simulaciones diferentes y comparamos el tiempo promedio entre eventos y la desviación estándar del tiempo entre eventos.
# Simular múltiples realizaciones con los parámetros estimados
simulaciones <- replicate(1000, {
sim <- simular_hawkes_exponencial(
N = 10,
lambda_inicial = lambda_0,
a = 0.2755 ,
beta = 1.3664 ,
b = 0.4958
)
return(sim$tiempos[-1])
})
# Comparar estadísticas
estadisticas_original <- c(
mean(diff(P)), # Tiempo promedio entre eventos observados
sd(diff(P)) # Desviación estándar observada
)
estadisticas_simuladas <- apply(simulaciones, 2, function(tiempos) {
c(mean(diff(tiempos)), sd(diff(tiempos)))
})
# Ver si las estadísticas son similares
texto_resumen <- sprintf("
================ RESULTADOS COMPARATIVOS ================
Estadísticas originales (datos observados):
- Tiempo medio entre eventos: %.4f
- Desviación estándar: %.4f
Estadísticas simuladas (promedio sobre %d simulaciones):
- Tiempo medio entre eventos: %.4f
- Desviación estándar: %.4f
=========================================================
",
estadisticas_original[1],
estadisticas_original[2],
ncol(simulaciones),
rowMeans(estadisticas_simuladas)[1],
rowMeans(estadisticas_simuladas)[2]
)
cat(texto_resumen)##
## ================ RESULTADOS COMPARATIVOS ================
##
## Estadísticas originales (datos observados):
## - Tiempo medio entre eventos: 2.4976
## - Desviación estándar: 2.8791
##
## Estadísticas simuladas (promedio sobre 1000 simulaciones):
## - Tiempo medio entre eventos: 2.3561
## - Desviación estándar: 2.7702
##
## =========================================================En efecto, vemos una considerable similitud entre las estadísticas calculadas, lo cual es un buen indicio del funcionamiento del método MCMC y de paso de la simulación del proceso de Hawkes.
También agregamos los tiempos dados por las 10 últimas simulaciones, para darle más profundidad a la comporación con los tiempos originales.
# Mostrar comparación entre datos reales y últimas 10 simulaciones
num_mostrar <- 10
ultimas_sim <- simulaciones[, (ncol(simulaciones) - num_mostrar + 1):ncol(simulaciones)]
lista_ultimas <- lapply(1:ncol(ultimas_sim), function(i) ultimas_sim[, i])
# Crear texto limpio para reporte
texto_comp <- "\n=========== COMPARACIÓN DE TIEMPOS OBSERVADOS Y SIMULADOS ===========\n\n"
# Agregar tiempos reales
texto_comp <- paste0(
texto_comp,
"Tiempos reales observados (P):\n ",
paste(round(P, 4), collapse = ", "),
"\n\nÚltimas ", num_mostrar, " simulaciones:\n"
)
# Agregar cada simulación
for(i in seq_along(lista_ultimas)) {
tiempos_i <- paste(round(lista_ultimas[[i]], 4), collapse = ", ")
texto_comp <- paste0(texto_comp,
sprintf("\n Simulación %d:\n %s\n", i, tiempos_i)
)
}
texto_comp <- paste0(texto_comp,
"\n======================================================================\n"
)
cat(texto_comp)##
## =========== COMPARACIÓN DE TIEMPOS OBSERVADOS Y SIMULADOS ===========
##
## Tiempos reales observados (P):
## 0.5532, 2.038, 3.861, 9.1488, 9.4211, 9.4344, 10.85, 10.8574, 19.4724, 23.0319
##
## Últimas 10 simulaciones:
##
## Simulación 1:
## 7.2508, 11.6568, 11.9023, 14.0822, 14.2316, 14.2847, 20.2287, 21.1804, 26.9402, 30.3905
##
## Simulación 2:
## 0.4942, 0.6862, 1.1712, 8.2485, 11.1177, 13.2856, 13.3211, 18.1565, 20.7989, 22.7237
##
## Simulación 3:
## 0.3892, 1.384, 1.4261, 1.4378, 1.4457, 3.4713, 5.6117, 6.1212, 7.9875, 9.869
##
## Simulación 4:
## 0.0446, 0.299, 0.4676, 0.6493, 1.0704, 1.7835, 14.259, 16.7971, 17.4182, 23.1772
##
## Simulación 5:
## 3.4864, 3.5398, 6.2095, 11.1118, 11.4788, 12.564, 15.1861, 15.5512, 20.2549, 23.6604
##
## Simulación 6:
## 0.1857, 0.2446, 1.3624, 9.3928, 11.5814, 12.8842, 13.1381, 16.71, 26.3836, 26.7437
##
## Simulación 7:
## 0.6623, 1.1647, 14.7598, 17.1882, 22.4796, 34.0194, 34.2837, 34.3627, 35.2365, 35.4961
##
## Simulación 8:
## 0.0972, 2.8215, 3.026, 3.5889, 3.7188, 4.5054, 4.7505, 7.9934, 11.3119, 12.8856
##
## Simulación 9:
## 0.6475, 0.6874, 0.9429, 1.3044, 2.2609, 3.932, 13.4412, 20.2842, 30.4563, 30.8705
##
## Simulación 10:
## 0.5002, 0.7766, 9.4067, 16.3107, 18.315, 19.3913, 26.8985, 38.3518, 39.8773, 51.3107
##
## ======================================================================Aquí podemos observar el fenómeno que expresamos antes. Si comparamos con una simulación específica como la “Simulación 10”, podría aparentar que el método MCMC fue un fracaso y la estimación es errónea, pero una vez tomamos la gran imagen de las simulaciones (con las estadísticas que hallamos) es que se entiende la efectividad de las estimaciones dadas.
Como mencionamos al inicio de esta sección de experimentación, podemos realizar este proceso en la dirección contraria para medir la efectividad del algoritmo. Es decir, podemos tomar parámetros fijos conocidos, simular un proceso de Hawkes con estos parámetros y luego con los tiempos obtenidos usar el algoritmo de Metropolis-Hastings para estimar de nuevo dichos parámetros. En la vida real, este proceso no teien mucho sentido, puesto que para qué buscaríamos estimar parámetros que ya conocemos, pero consideramos que es un buen experimento para de nuevo medir la efectividad del método.
Siguiendo el ejemplo anterior, tomamos los párametros de acuerdo a los que estimamos en el ejemplo anterior y hacemos repetidas veces su respectiva simulación. A partir de estas simualciones, tomamos unos tiempos representativos lo cuales usaremos luego para aplicar el método MCMC:
# Parámetros
a <- 0.2755
beta <- 1.3664
b <- 0.4958
lambda_0 <- 1.5
T_max <- 100
# Número de simulaciones
R <- 100000
# Lista para almacenar todas las simulaciones
todas_las_simulaciones <- list()
for (r in seq_len(R)) {
sim <- simular_hawkes_exponencial(
N = 10,
lambda_inicial = lambda_0,
a = a,
beta = beta,
b = b
)
# Extraer tiempos menores a T_max
t <- sim$tiempos
t <- t[t <= T_max]
# Quitar T=0 si aparece
if (length(t) > 1 && t[1] == 0) t <- t[-1]
todas_las_simulaciones[[r]] <- t
}
# Encontrar simulación "típica"
num_eventos <- sapply(todas_las_simulaciones, length)
num_events_typical <- round(mean(num_eventos))
sim_representativa <- todas_las_simulaciones[[which.min(abs(num_eventos - num_events_typical))[1]]]
tiempos_representativos <- sim_representativa
# Salida limpia en una sola impresión
salida <- paste0(
"Número de eventos en simulación representativa: ", length(tiempos_representativos), "\n",
"Primeros 10 tiempos: ",
paste(head(tiempos_representativos, 10), collapse = ", ")
)
cat(salida, "\n")## Número de eventos en simulación representativa: 10
## Primeros 10 tiempos: 0.131548378799696, 1.09993420464854, 2.03791354979777, 7.23536754963583, 17.1838304393611, 18.2108683661392, 18.3370583268058, 18.7539838236615, 21.6520845570306, 27.9065862181863Ahora, usando estos tiempos representativos, y con valores iniciales apropiados, aplicamos el método MCMC:
# Tiempos representativos
c(0.0444679305014493,
0.325817576948851,
5.36830635607417,
14.4773172574539,
18.1374117640713,
18.5499058218557,
25.7638200865306,
26.7101940708983,
31.5536866363207,
34.9438796047865)## [1] 0.04446793 0.32581758 5.36830636 14.47731726 18.13741176 18.54990582
## [7] 25.76382009 26.71019407 31.55368664 34.94387960# Parámetros razonables
lambda_0 <- 1.5
a_init <- 0.8
b_init <- 0.5
beta_init <- 1.2
T_final <- max(P) + 5
n_iter <- 100000
burn_in <- 10000
# Sigmas ajustados
sigma_a <- 0.3
sigma_b <- 0.3
sigma_beta <- 0.7
salida <- capture.output(
res <- metropolis_hastings_hawkes(
P = P, n_iter = n_iter, burn_in = burn_in,
a_init = a_init, b_init = b_init, beta_init = beta_init,
sigma_a = sigma_a, sigma_b = sigma_b, sigma_beta = sigma_beta,
T_final = T_final, lambda_0 = lambda_0
)
)
salida_filtrada <- salida[!grepl("Iteración", salida)]
cat(salida_filtrada, sep = "\n")## Ejecutando MCMC...
## RESULTADOS MCMC - PROCESO DE HAWKES
## ESTIMACIONES:
## a (tasa base): 0.2634
## b (escala saltos): 0.5372
## beta (decaimiento): 1.3753
##
## TASAS DE ACEPTACIÓN:
## a: 0.382
## b: 0.451
## beta: 0.509De nuevo, podemos apreciar que la estimación dada por el algoritmo de Metropolis-Hastings es satisfactoria, sin ser idéntica. Pensamos que esto puede suceder debido a que no es sencillo dar unos tiempos que representen de manera ideal el proceso que se quiere simular, y que las estimaciones no sean idénticas a los parámetros reales no es sinónimo de que estas no puedan generar un proceso de Hawkes con tiempos cercanos a los observados, como vimos en la experimentación anterior
Se concluye que el método es una opción viable en cualquiera de las dos situaciones que ilustramos en esta sección de experimentación.
[1] Laub, P. J., Lee, Y., & Taimre, T. (2021). The Elements of Hawkes Processes. Springer International Publishing. https://doi.org/10.1007/978-3-030-84639-8
[2] Casella, G., & Robert, C. P. (2004). Monte Carlo Statistical Methods (2nd ed.). Springer-Verlag. https://doi.org/10.1007/978-1-4757-4145-2
[3] Dassios, A., & Zhao, H. (2013). Exact simulation of Hawkes process with exponentially decaying intensity. Electronic Communications in Probability, 18(62), 1-13. https://doi.org/10.1214/ECP.v18-2717
[4] Ozaki, T. (1979). Maximum likelihood estimation of Hawkes self-exciting point processes. Annals of the Institute of Statistical Mathematics, 31(1), 145–155. DOI: 10.1007/BF02480272