Simulación Exacta por Composición

Se parte del hecho que se quiere simular un proceso de Hawkes \(\{T_i\}_{i\in\mathbb{N}}\), con un kernel de decaimiento exponencial, donde su función de intensidad condicional se expresa de la siguiente forma:

\[ \lambda^*(t)=\lambda+\sum_{t_i<t} \alpha e^{-\beta(t-t_i)} \]

Donde los parámetros son \(\lambda\) la intensidad base, \(\alpha\) coeficiente de auto-excitación y \(\beta\) tasa de decaimiento exponencial.

Los autores Dassios y Zhao en 2013 [1], proponen un algoritmo de simulación exacta, usando la propiedad de Márkov del kernel exponencial. Para esto se reescribe la función de intensidad condicional, de forma que depende solo del evento que sucede antes de \(t\), tal que \(t_k<t<t_{k+1}\), es decir \(t_k\) es el último evento antes de \(t\).

\[ \lambda^*(t)=\lambda+(\lambda^*(t_k)+\alpha-\lambda)e^{\beta t_k}e^{-\beta t} \]

En el artículo, se llega al resultado de que a partir de dos variables \(T^{(1)}_{k+1}\) y \(T^{(2)}_{k+1}\), se puede obtener la distribución de \(T_{k+1}\), a través del método de la transformada inversa. Los resultados de presentan a continuación.

Resultado 1: Si \(T^{(1)}_{k+1}\) y \(T^{(2)}_{k+1}\), son variables aleatorias independientes y poseen las siguientes distribuciones.

\[\begin{align*} P(T_{k+1}^{(1)} > t) &= \exp\left( -(t - t_k)\lambda \right) \\ P(T_{k+1}^{(2)} > t) &= \exp\left( - (\lambda^*(t_k) + \alpha - \lambda)\beta^{-1} \left[1 - e^{-\beta(t-t_k)} \right] \right) \end{align*}\]

Entonces, se cumple que: \[ min\{T_{k+1}^{(1)},T_{k+1}^{(2)}\}=T_{k+1} \]

Donde \(T_{k+1}\) es un proceso de Hawkes de parámetros \(\lambda\) , \(\alpha\) y \(\beta\).

A partir de este resultado, se plantea un algoritmo para obtener simulaciones de un proceso de Hawkes, los pasos del algoritmo se detallan a continuación:

Algoritmo

  1. Fije sus parámetros \(\lambda\) , \(\alpha\) y \(\beta\) y escoja \(N\) el número de observaciones que quiere simular

  2. Inicializar el algoritmo con: \[\begin{align*} t_0 &= 0 \\ \lambda^*(t_0) &= \lambda^*(0) = \lambda \end{align*}\]

Inicie un ciclo for de forma que el tiempo \(k+1\) se obtiene de la siguiente forma

  1. Obtenga muestras aleatorias independiente de dos uniformes, tal que:

\[\begin{align*} U^{(1)} &\sim U(0,1) \\ U^{(2)} &\sim U(0,1) \end{align*}\]

  1. Evalué en la inversa de las distribuciones acumuladas de \(T^{(1)}_{k+1}\) y \(T^{(2)}_{k+1}\), de la siguiente forma que:

\[\begin{align*} T_{k+1}^{(1)} &= t_k - \lambda^{-1} \ln(U^{(1)}) \\ T_{k+1}^{(2)} &= t_k - \beta^{-1} \ln \left( 1 + \beta (\lambda^*(t_k) + \alpha - \lambda)^{-1} \ln(U^{(2)}) \right) \end{align*}\]

  1. Tome el mínimo entre las dos variables aleatorias de forma que:

\[ min\{T_{k+1}^{(1)},T_{k+1}^{(2)}\}=T_{k+1} \]

  1. Almacene el tiempo en que ocurre el evento \(k+1\) y calcule y almacene \(\lambda^*(t_{k+1})\), para la siguiente iteración del ciclo.

Fin del ciclo for

  1. Repita hasta obtener \(N\) observaciones del proceso.

Implementación en R

Se propone la siguiente función que ejecuta, los pasos anteriormente mencionados.

SimHawkesExact = function(alpha, beta, lambda, N){
  #Vector Donde Guardar Procesos Puntual
  T_times     <- numeric(N)
  
  #Inicializar
  lambda_tk   <- lambda
  tk          <- 0
  
  #Algoritmo
  
  for (i in 1:N) {
    #Generar dos numeros aleatorios entre 0 y 1
    
    U1 <- runif(1)
    U2 <- runif(1)
    
    # T^(1)_{k+1}
    T1 <- tk - (1 / lambda) * log(U1) 
    
    # T^(2)_{k+1}
    A <- lambda_tk + alpha - lambda 
    arg <- 1 + (beta / A) * log(U2)

    if (arg<=0 ){
      T2=Inf}
      else{
        T2 <- tk - (1 / beta) * log(arg)
      }
    
    
    # siguiente tiempo y actualización de intensidad según 
    
   
    t_next <- min(T1, T2)
    
    T_times[i] <- t_next
    
    #Actualizar los parametros para la siguiente iteración
    lambda_tk= lambda+(lambda_tk+alpha-lambda)*exp(-beta*(t_next-tk))
    tk= t_next
    
    
    
    
  }
return(T_times)
  
}

Como ejemplo se van a generar 50 observaciones de un proceso de Hawkes con párametros \(\lambda=7\) , \(\alpha=1.5\) y \(\beta=3\), recordemos que se debe cumplir que \(\alpha<\beta\) para que el proceso no explote.

set.seed(123)
alpha = 1.5
beta=3
lambda=7
Simulacion=SimHawkesExact(alpha = alpha, beta=beta, lambda=lambda,N=50)
head(Simulacion)
## [1] 0.1780375 0.2349144 0.2436827 0.2676839 0.3527169 0.3590207

Adicionalmente, se grafica la función de intensidad condicional, para observar cómo se comporta el proceso.

#Graficar la intensidad
# Función de intensidad λ*(t) dada la historia de eventos
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))))
  }
}


t <- seq(0, max(Simulacion) + 0.5, length.out = 5000)


lambda_vals <- sapply(t, function(x) inten(x, eventos=Simulacion, alpha=alpha, beta=beta, lambda=lambda))


plot(t, lambda_vals, type = "l", lwd = 1.5, col = "#6495ED", # Un azul bonito
     xlab = "t", ylab = expression(paste(lambda, "*(t)")),,
     main = "Intensidad Condicional de un Proceso de Hawkes",
     ylim = c(0, max(lambda_vals)*1.1)) # Un poco de aire arriba



abline(h = lambda, col="gray", lty=2) # Línea de la intensidad base
legend("topright", legend=c("Intensidad Base"), 
       col=c("gray"), lty=c(2), cex=0.8)

Estimació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)
set.seed(12)
Simulacion=SimHawkesExact(alpha = 11.3, beta=44.1, lambda=22.7,N=10000)

# 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 
mu_init    <- 15         
alpha_init <- matrix(7, nrow = 1, ncol = 1)
beta_init  <- matrix(30, 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 = "Nelder-Mead"
)

## 5. Ver resumen y parámetros estimados
theta=fit$estimate
summary(fit)  
## --------------------------------------------
## Maximum Likelihood estimation
## BFGS maximization, 23 iterations
## Return code 0: successful convergence 
## Log-Likelihood: 24414.43 
## 3  free parameters
## Estimates:
##        Estimate Std. error t value Pr(> t)    
## mu1    22.34069    0.08343   267.8  <2e-16 ***
## alpha1 12.04933    0.06751   178.5  <2e-16 ***
## beta1  46.04821    0.06514   706.9  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------

Nota: 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)
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  Delta
## D = 0.0054369, p-value = 0.929
## 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 = 0.21992, p-value = 0.984
## alternative hypothesis: NA

El 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)
lag.plot(U,main="U_(k-1) vs U_k")

acf(U, lag.max = 10, main = "ACF de U")

Box.test(U, lag = 20, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  U
## X-squared = 30.537, df = 20, p-value = 0.06161

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

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

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

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

  2. 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 \]

  1. 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)
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  S
## D = 0.0057178, p-value = 0.8993
## alternative hypothesis: two-sided

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

  1. 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].

  2. Aplicar la transformación \(M(t)\), a la secuencia construida anteriormente y hallar el punto en el tiempo donde alcanza su máximo \(M^*\).

  3. 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.556346171 0.998458667

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

Prueba del tiempo 1

Se sabe que si \(M(t)\) es un movimiento Browniano, entonces \(M(1)\sim N(0,1)\), entonces la prueba consta en ver si \(M(1)\) cae en el intervalo \([Z_{\alpha /2},Z_{1-\alpha /2}]\).

Implementación en R

Esta prueba, no está en ningún paquete, por lo que se implementa de manera manual.

BMTHawkesSimple = function(RTCT,sig=0.05){
  t = RTCT/max(RTCT)
  nt=seq(1,length(RTCT))
  Mt=(nt-t*max(RTCT))/sqrt(max(RTCT))
  M1=Mt[length(RTCT)]
  Li=qnorm(p=sig/2)
  Ls=qnorm(p=1-(sig/2))
  return(c(Li,M1,Ls))
  
}
BMTHawkesSimple(RTCT)
## [1] -1.959964e+00 -5.789507e-05  1.959964e+00

Se observa que \(M(1)\) 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.

Referencias

[1] Dassios, A., & Zhao, H. (2013). Exact simulation of Hawkes process with exponentially decaying intensity. Electronic Communications in Probability, 18(62), 1–13.

[2] Ozaki, T. (1979). Maximum likelihood estimation of Hawkes’ self-exciting point processes. Annals of the Institute of Statistical Mathematics, 31(1), 145–155.

[3] Ogata, Y. (1978). The asymptotic behaviour of maximum likelihood estimators for stationary point processes. Annals of the Institute of Statistical Mathematics, 30(1), 243–261.

[4] Laub, P. J., Lee, Y., & Taimre, T. (2021). The Elements of Hawkes Processes. Springer Cham. https://doi.org/10.1007/978-3-030-84639-8

[5] Lee, K. (2025). emhawkes: Exponential Multivariate Hawkes Model. R package version 0.9.8. https://CRAN.R-project.org/package=emhawkes

[6] Lewis, P. A. W. (1965). Some results on tests for Poisson processes. Biometrika, 52(1/2), 67–77.

[7] Kim, S.-H., & Whitt, W. (2015). The Power of Alternative Kolmogorov-Smirnov Tests Based on Transformations of the Data. ACM Transactions on Modeling and Computer Simulation (TOMACS), 25(4), Article 24.