Procesos de Hawkes en Acciones de Ecopetrol
Sebastian Romero & Julian Marulanda
2025-12-07
Carga de Datos
datos=read.csv("C:\\Users\\User\\Desktop\\Procesos Actuaria\\ecopetrol_hoy_hawkes.csv")
head(datos)## Fecha_Original T_minutos Volumen
## 1 2025-12-05 14:30:00+00:00 0.773956 18300
## 2 2025-12-05 14:31:00+00:00 1.438878 786
## 3 2025-12-05 14:32:00+00:00 2.858598 7293
## 4 2025-12-05 14:33:00+00:00 3.697368 7747
## 5 2025-12-05 14:34:00+00:00 4.094177 5451
## 6 2025-12-05 14:35:00+00:00 5.975622 5839summary(datos$Volumen)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 100 1194 2551 7492 7764 114492fecha_inicio <- datos$Fecha_Original[1]
fecha_fin <- datos$Fecha_Original[nrow(datos)]
T_maximo <- tail(datos$T, 1)
datosmayores=subset(datos, datos$Volumen>=7764)
head(datosmayores) ## Fecha_Original T_minutos Volumen
## 1 2025-12-05 14:30:00+00:00 0.773956 18300
## 8 2025-12-05 14:37:00+00:00 7.786064 12788
## 12 2025-12-05 14:41:00+00:00 11.926765 16828
## 57 2025-12-05 15:30:00+00:00 60.458916 20238
## 74 2025-12-05 15:48:00+00:00 78.233939 25723
## 80 2025-12-05 15:54:00+00:00 84.783898 23966Simulacion=datosmayores$T_minutos
summary(datos$Volumen)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 100 1194 2551 7492 7764 114492Estimación de Páramaetros por Máxima Verosimilitud
Ozaki en 1979 [2], dedujo matemáticamente formas para la log verosimilitud de un proceso de Hawkes con decaimiento exponencial, de igual forma obtuvo formas para el gradiente y la matriz Hessiana, de tal forma que se pueden aplicar algoritmos de la forma Newton-Rapshon. Los resultados fueron los siguientes:
Resultado 2: Se observa un proceso de Hawkes hasta el tiempo \(T\), de forma que se observan \(n(T)\) eventos, entonces se tienen unas observaciones de la forma \(\{t_1,t_2,t_3,...,t_{n(T)}\}\). Entonces su log verosimilitud es de la forma.
\[ A(i) = e^{-\beta(t_i - t_{i-1})}(1 + A(i-1)) \hspace{0.7cm} A(1)=0 \] \[ \ell = \sum_{i=1}^{n(T)} \ln(\lambda + \alpha A(i)) - \lambda T - \frac{\alpha}{\beta} \sum_{i=1}^{n(T)} \left( 1 - e^{-\beta(T-t_i)} \right) \]
Adicionalmente, Ogata en 1978 [3] demostró las propiedades asintóticas de la estimaciones máximo verosímiles de un proceso puntual estacionario. La estacionariedad de un proceso de Hawkes está demostrada en Laub, Lee y Taimre 2021 [4]. Cabe aclarar que se hace referencia a un proceso estacionario débil.
Resultado 3: Si \(\hat{\boldsymbol{\theta}}_T= (\hat\lambda , \hat\alpha, \hat\beta)\), son estimadores máximo verosímiles, para un proceso de Hawkes con decaimiento exponencial, entonces cumplen las siguientes propiedades.
Convergencia Casi Siempre: El estimador converge casi siempre al verdadero valor de los parámetros: \[ \hat{\boldsymbol{\theta}}_T \xrightarrow{c.s.} \boldsymbol{\theta} \]
Normalidad Asintótica: La distribución del estimador cuando \(T \to \infty\): \[ \sqrt{T} \left( \hat{\boldsymbol{\theta}}_T - \boldsymbol{\theta} \right) \xrightarrow{d} \mathcal{N} \left( \mathbf{0}, \mathcal{I}^{-1}(\boldsymbol{\theta}) \right) \] O equivalentemente: \[ \hat{\boldsymbol{\theta}}_T \xrightarrow{d} \mathcal{N} \left( \boldsymbol{\theta}, \frac{\mathcal{I}^{-1}(\boldsymbol{\theta})}{T} \right) \] Donde \(\mathcal{I}(\boldsymbol{\theta})\) es la Matriz de Información de Fisher
Eficiencia Asintótica: El estimador es asintóticamente eficiente, lo que significa que su varianza asintótica alcanza la cota inferior de Cramér-Rao.
Implementación en R
El paquete emhawkes [5] ya implementa este método, entonces se hace una simulación con el algoritmo de composición y se estiman los parámetros, se usan valores de \(\alpha=11.3\) , \(\beta=44.1\) y \(\lambda=22.7\)
library(emhawkes)
# inter-arrival: primer valor 0, luego las diferencias entre tiempos
inter_arrival <- c(0, diff(Simulacion))
# todo pertenece a la misma dimensión (univariado)
type <- rep(1L, length(Simulacion))
## 3. Definir el modelo
# Valores iniciales para los parámetros
intensidad_promedio <- length(Simulacion) / max(Simulacion)
# 2. Usamos esta información para los valores iniciales
# Asumimos que la mitad es base (mu) y la mitad es excitación (Hawkes)
mu_init <- intensidad_promedio / 2
alpha_init <- matrix(0.5, nrow = 1, ncol = 1)
beta_init <- matrix(1.0, nrow = 1, ncol = 1)
# Crear el objeto hspec (modelo)
h <- new("hspec",
mu = mu_init, # baseline
alpha = alpha_init, # matriz 1x1
beta = beta_init) # matriz 1x1
## 4. Ajustar el modelo por máxima verosimilitud con hfit
fit <- hfit(
object = h,
inter_arrival = inter_arrival,
type = type,
method = "NM", # <--- CAMBIO AQUÍ: Nelder-Mead
control = list(iterlim = 5000) # Damos más iteraciones
)
## 5. Ver resumen y parámetros estimados
summary(fit) ## --------------------------------------------
## Maximum Likelihood estimation
## Nelder-Mead maximization, 174 iterations
## Return code 0: successful convergence
## Log-Likelihood: -200.2841
## 3 free parameters
## Estimates:
## Estimate Std. error t value Pr(> t)
## mu1 0.05164 0.02337 2.209 0.027144 *
## alpha1 0.08172 0.02341 3.491 0.000482 ***
## beta1 0.10001 0.02745 3.643 0.000270 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------#Lambda, alpha y beta orden
theta=fit$estimateNota: El paquete toma \(\mu=\lambda\)
Bondad de Ajuste del Modelo
En el libro de Laub, Lee y Taimre 2021 [4]. Se dedica una sección entera a la bondad de ajuste, acá se muestran esos metodos y su implementación en R.
Varias de estas pruebas de basan en los siguientes dos teoremas:
Teorema 1 (Teorema del Cambio de Tiempo Aleatorio)
Sea \(\{t_1, t_2, \dots, t_k\}\) una realización en el intervalo de tiempo \([0, T]\) de un proceso puntual con función de intensidad condicional \(\lambda^*(\cdot)\).
Si \(\lambda^*(\cdot)\) es positiva sobre \([0, T]\) y \(\Lambda(T) < \infty\) casi seguramente, entonces los puntos transformados \(\{\Lambda(t_1), \Lambda(t_2), \dots, \Lambda(t_k)\}\) forman un proceso de Poisson con tasa unitaria.
Teorema 2 (Análisis del Residuo)
Considere una secuencia no acotada, creciente monótona de puntos en el tiempo \(\{t_1,t_2....\}\) y un compensador monótono, continuo \(\Lambda(\cdot)\) tal que \(lim_{t \to \infty}\Lambda(t)= \infty\) casi seguramente.
Entonces la secuencia transformada \(\{t^*_1,t^*_2....\}=\{\Lambda(t_1),\Lambda(t_2)....\}\) que proviene del proceso de conteo \(N^*(t)\), es un proceso de Poisson con tasa 1 si y solo si la secuencia \(\{t_1,t_2....\}\) proviene de un proceso puntual definido por \(\Lambda(\cdot)\).
Compensador de un procesos de Hawkes con kernel exponencial
\[ \Lambda(t) = \int_{0}^{t} \lambda^{*}(s)\, ds = \lambda t + \frac{\alpha}{\beta} \sum_{t_i < t} \left[ 1 - e^{-\beta(t - t_i)} \right] \]
Pruebas Básicas
Se calculan los tiempos entre llegadas de la secuencia transformada por el compensador, es decir:
\[ \{\tau_1,\tau_2,\tau_3...\}=\{t^*_1,t^*_2-t^*_1,t^*_3-t^*_2....\} \]
Se hace un qq-plot, contra los cuantiles de una exponencial con tasa 1 o también se pueden hacer pruebas como Kolmogrov-Smirnov o Anderson-Darling.
Implementación en R
hawkes_compensator_exp <- function(t, H_t, lambda, alpha, beta) {
if (t <= 0) return(0)
# Parte determinista: λ * t
comp <- lambda * t
# Asegurarnos de usar solo eventos estrictamente antes de t (opcional pero sano)
H_t <- H_t[H_t < t]
if (length(H_t) > 0) {
comp <- comp + sum((alpha / beta) * (1 - exp(-beta * (t - H_t))))
}
return(comp)
}
RTCT=numeric(0)
for (i in 1:length(Simulacion)) {
RTCT[i] = hawkes_compensator_exp(Simulacion[i],Simulacion,theta[1],theta[2],theta[3])
}
Delta <- diff(c(0, RTCT))
#Guardar los parametros ajustados
#install.packages("EnvStats")
library(EnvStats)##
## Adjuntando el paquete: 'EnvStats'## The following objects are masked from 'package:stats':
##
## predict, predict.lm## The following object is masked from 'package:base':
##
## print.default# QQ-plot directamente
EnvStats::qqPlot(Delta, dist = "exp",
param.list = list(rate = 1),
add.line = TRUE)
#Kolmogroc smirnov
ks.test(Delta, "pexp", rate = 1)##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: Delta
## D = 0.11914, p-value = 0.1353
## alternative hypothesis: two-sided#Anderson Darling
library(ADGofTest)
ad.test(Delta, distr.fun = pexp, rate = 1)##
## Anderson-Darling GoF Test
##
## data: Delta and pexp
## AD = 1.9912, p-value = 0.093
## alternative hypothesis: NAEl qqplot, muestra un comportamiento muy bueno y los test no rechazan \(H_0\) es decir no se rechaza que provengan de una distribución exponencial con tasa 1.
Pruebas de Independencia
Los \(\tau_i\) deben ser independientes uno del otro, entonces se pueden hacer la prueba de Ljung-Box para ver la significancia de las correlaciones conjuntamente, prueba de Barlett para ver la significancia de cada correlación. el grafico de rezago y autocorrelograma de las variables \(U_i\), que se definen como:
\[ U_k=1-e^{-(t^*_k-t^*_{k-1})} \]
Nota:El autocorrelograma nos permite hacer la prueba de Bartlett para todos los rezagos de una manera grafica con la línea punteada azul.
Implementación en R
U <- 1 - exp(-Delta)
stats::lag.plot(U,main="U_(k-1) vs U_k",labels = F,do.lines = F)
acf(U, lag.max = 10, main = "ACF de U")
Box.test(U, lag = 20, type = "Ljung-Box")##
## Box-Ljung test
##
## data: U
## X-squared = 23.54, df = 20, p-value = 0.263No se rechaza \(H_0\) de la prueba de Ljung-Box es decir no se rechaza que todas las correlaciones hasta el rezago 20 sean 0 y en las pruebas graficas tampoco se observan patrones o correlaciones altas.
Prueba de Lewis
Propuesta por Lewis en 1965 [6] y descrita en detalle en Kim y Whitt en 2015 [7], es una prueba con mayor potencia, que la de Kolmogrov-Smirnov, en los casos donde se quiere probar si unos datos provienen de una variable exponencial con tasa 1. El método de la prueba se describe a continuación.
- Se calcula una nueva secuencia de observaciones de la siguiente forma:
\[ \{t^*_1/t^*_N,t^*_2/t^*_N,...,t^*_{N-1}/t^*_N\}=\{U_{(1)},U_{(2)},...,U_{(n)}\} \]
- Se calculan las variables aleatorias \(C_i\) de la siguiente forma:
\[ C_1=U_{(1)} \hspace{1cm} C_j=U_{(j)}-U_{(j-1)} \hspace{1cm} C_{n+1}=1-U_{(n)} \]
Ahora se ordenan las variables \(C_i\) de forma que \(C_{(1)}\) es el valor más pequeño y se cumple que: \[ C_{(1)}<C_{(2)}<...<C_{(n+1)} \]
Ahora se construyen nuevas variables \(Z_j\) a partir de las anteriores de forma que:
\[ Z_j=(n+2-j)(C_{(j)}-C_{(j-1)}) \hspace{1cm} C_{(0)}=0 \]
- Por último se calcula las sumas acumuladas de \(Z_j\) y se hace el test Kolmogrov-Smirnov sobre las variables \(S_k\), comparándolas con una \(U(0,1)\).
\[ S_k=Z_1+Z_2+...+Z_k \hspace{0.5cm} k\in\{1,2,3,...,n\} \]
Implementación en R
Esta prueba, no está en ningún paquete, por lo que se implementa de manera manual.
LewisTest= function(RTCT){
U = RTCT/max(RTCT)
U = U[-length(U)]
# n
n = length(U)
# 2. Calcular Intervalos C
C_raw = diff(c(0, U, 1))
#3 Ordenar las C
C = sort(C_raw)
# 4. Calcular las diferencias de los intervalos ordenados
diffs_C = c(C[1], diff(C))
# 5. Aplicar la ponderación
Z = numeric(length(diffs_C))
for (j in 1:length(Z)) {
Z[j] = diffs_C[j] * (n + 2 - j)
}
#6 Suma Acumulada
S = cumsum(Z)
S=S[1:n]
resultado = ks.test(S, "punif", min=0, max=1)
return(resultado)
}
LewisTest(RTCT)##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: S
## D = 0.16177, p-value = 0.01508
## alternative hypothesis: two-sidedNo se rechaza \(H_0\), por lo que no se rechaza que los datos provengan de una exponencial con tasa 1.
Pruebas Basadas en el Movimiento Browniano
Todas estas pruebas y sus descripciones se encuentran en Laub, Lee y Taimre 2021 [4].
Se tiene que \(n(t)\) es un proceso de Poisson con tasa \(T\) y se define \(M(t)\) como: \[ M(t)=\frac{n(t)-tT}{\sqrt{T}} \] Por el principio de invarianza de Donsker se tiene que si \(T \to \infty\), entonces \((M(t) : t\in[0,1])\) converge en distribución a un movimiento Browniano estandar \((B(t) : t\in[0,1])\).
Pruebas de la Ley del Arcoseno
La ley del arcoseno dice que el punto en el tiempo \(M^*\) donde el movimiento Browniano alcanza su máximo en el intervalo [0,1], distribuye Beta(1/2,1/2). Formalmente se nota como:
\[ M^* = \operatorname*{arg\,max}_{s \in [0,1]} B(s) \sim Beta(1/2,1/2) \]
Con este resultado los autores proponen la siguiente prueba.
Construir la secuencia \(\{t^*_1/\Lambda (T),t^*_2/\Lambda (T),...,t^*_k/\Lambda (T)\}\), donde \(t_k\) es el ultimo evento observado y \(T\) al tiempo hasta que fue observado el proceso, esta secuencia debería seguir una distribución Poisson con tasa \(T\), sobre el intervalo [0,1].
Aplicar la transformación \(M(t)\), a la secuencia construida anteriormente y hallar el punto en el tiempo donde alcanza su máximo \(M^*\).
No rechazar con una significancia de \(\alpha\) que las secuencia proviene de un proceso de Poisson con tasa 1, si \(M^*\) cae el intervalo generado por los cuantiles \((\alpha/2,1-\alpha/2)\) de una distribución \(Beta(1/2,1/2)\).
Implementación en R
Esta prueba, no está en ningún paquete, por lo que se implementa de manera manual.
BMTHawkes = function(RTCT,sig=0.05){
t = RTCT/max(RTCT)
nt=seq(1,length(RTCT))
Mt=(nt-t*max(RTCT))/sqrt(max(RTCT))
idx_max = which.max(Mt)
Mest = t[idx_max]
Li=qbeta(p=sig/2,shape1=1/2,shape2=1/2)
Ls=qbeta(p=1-(sig/2),shape1=1/2,shape2=1/2)
return(c(Li,Mest,Ls))
}
BMTHawkes(RTCT)## [1] 0.001541333 0.404352587 0.998458667Se observa que \(M^*\) cae en el intervalo, por lo tanto, no se rechaza \(H_0\) es decir no se rechaza que la secuencia proviene de un proceso de Poisson con tasa 1.
# 1. Definir la función de intensidad (Igual que antes)
inten <- function(t, eventos, alpha, beta, lambda) {
eventos_pasados <- eventos[eventos < t]
if (length(eventos_pasados) == 0) {
return(lambda)
} else {
return(lambda + sum(alpha * exp(-beta * (t - eventos_pasados))))
}
}
# 2. Generar la secuencia de tiempo 't' en MINUTOS (para el cálculo matemático)
t <- seq(0, max(Simulacion) + 0.5, length.out = 1000)
# 3. Calcular la intensidad (Usando los minutos, porque los parámetros dependen de ello)
# Asegúrate de que theta tenga los valores estimados: theta = fit$estimate
lambda_vals <- sapply(t, function(x) inten(x, eventos=Simulacion, alpha = theta[2], beta = theta[3], lambda = theta[1]))
# --- PARTE NUEVA: CONVERSIÓN A HORAS PARA EL GRÁFICO ---
# A. Definir la hora de inicio (T=0).
# Pon la fecha de tus datos. La hora debe ser 14:30:00 (Apertura NYSE)
fecha_base <- as.POSIXct("2025-12-05 14:30:00", tz = "UTC")
# B. Crear el vector de tiempo real
# Sumamos t * 60 porque as.POSIXct suma segundos, y t está en minutos.
tiempo_real <- fecha_base + (t * 60)
# C. Graficar usando el tiempo real
plot(tiempo_real, lambda_vals, type = "l", lwd = 2, col = "red",
xlab = "Hora (UTC)",
ylab = expression(paste( lambda, "*(t)")),
main = "Intensidad Condicional Estimada",
# xaxt = "n" # Descomenta esto si quieres personalizar más el eje
)
# Opcional: Añadir una cuadrícula de tiempo para que se lea mejor
axis.POSIXct(1, x = tiempo_real, format = "%H:%M") # Formato Hora:Minuto
grid()