El algoritmo presentado en este documento tiene como objetivo simular un Proceso de Hawkes con tiempo entre eventos “espontáneos” (immigrants) exponencial con intensidad \(\lambda_0(t)\) y con kernel exponencial de la forma \(\phi(t)=\gamma e^{-\delta t}\) [1].

En particular, tomamos un proceso de Hawkes con intensidad

\[\lambda(t)=a+(\lambda_0-a) e^{-\delta t}+\sum_{T_i <t}\gamma e^{\delta (t-T_i)} \ , \ t>0\]

Las siguientes dos secciones son una síntesis e interpretación del trabajo hecho en [2].

Propiedad Markoviana

La función de intensidad del proceso de Hawkes presentado coincide con nuestra definición de proceso estocástico, considerando \(\{\lambda(t):t\geq0\}\) con conjunto índice \(T=[0,\infty)\). Además, debido a la elección del núcleo exponencial, cumple con la propiedad markoviana; esta será la base del algoritmo presentado más adelante. En efecto, es fácil comprobar que \(\lambda(t)\) satisface la EDO:

\[\frac{d\lambda(t)}{dt}=-\delta (\lambda(t)-a)\].

Tomando \(T_k \leq t <T_{k+1}\) (tiempo dado entre el \(k-\)ésimo y el \((k+1)\)-esimo evento del proceso) tendremos que la solución de la EDO es:

\[\lambda(t)=(\lambda(T_k^+)-a)e^{-\delta (t-T_k)}+a\].

Nota: Se toma la condición inicial \(\lambda(T_k^+)=\lambda(T_k)+\gamma\) (Salto inmediato en el evento).

Esto es justamente la propiedad de Markov, puesto que la intensidad del proceso inmediatamente después del último evento solo depende de su valor en el último evento.

Descomposición del proceso

Sea \(\tau\) el tiempo entre el \(k-\)ésimo evento y el \(k+1-\)ésimo evento, es decir, \(\tau=T_{k+1}-T_k\).

\[\begin{aligned} F_{\tau}(s) = \mathbb{P}(\tau \leq s) &= 1 - \mathbb{P}(\tau > s) \\ &= 1 - \mathbb{P}(N_{T+s} - N_T = 0) \quad \text{(Proceso de conteo)} \\ &= 1 - \exp\left(-\int_T^{T+s} \lambda(u)du\right) \quad \text{(Proceso de Poisson no homogéneo)} \\ &= 1 - \exp\left(-\int_0^{s} \lambda(T + v)dv\right) \quad \text{(Cambio de variable)} \\ &= 1 - \exp\left(-\int_0^{s} \left[a + (\lambda(T^+) - a)e^{-\delta v}\right] dv\right) \quad \text{(Propiedad de Markov)} \\ &= 1 - \exp\left(-as - (\lambda(T^+) - a)\frac{1 - e^{-\delta s}}{\delta}\right) \quad \text{(Evaluando la integral)} \end{aligned}\]

De esta manera, tenemos que \[\mathbb{P}(\tau>s)=\mathbb{P}(s_0>s) \mathbb{P}(s_1>s)\] con:

Esta es la descomposición de los tiempos entre llegadas, la cual usamos para muestrear valores de \(\tau\).

El siguiente algoritmo lo tomamos de [1]:

Algoritmo

  1. Tome \(T_0=0\) con \(\lambda(T_0^+)=\lambda(T_0^-)=\lambda_0\) tasa inicial.
  2. Para \(i=1,2 \ldots,N\)
    1. Tome \(u_0 \sim U(0,1)\) y defina \(s_0=-\frac{1}{a}\text{ln}(u_0)\).
    2. Tome \(u_1 \sim U(0,1)\) y defina \(d=1+\frac{\delta\text{ln}(u_1)}{\lambda(T_{i-1}^+)-a}\).
    3. Si \(d>0\), tome \(s_1=-\frac{1}{\delta}\text{ln}(d)\) y \(\tau_i=\text{mín}(s_0,s_1)\), e.o.c tome \(\tau_i=s_0\).
    4. Tome \(T_i=T_{i-1}+\tau_i\).
    5. Tome \(\lambda(T_i^-)=(\lambda(T_{i-1}^+)-a)e^{-\delta \tau_i}+a\).
    6. Tome \(\lambda(T_i^+)=\lambda(T_i^-)+\gamma\)

Implementación

Paso 1

Definimos los parámetros para nuestro proceso de Hawkes.

lambda_inicial <- 0.2   
alpha <- 0.6       
delta <- 1.2      
a <- 0.2
N <- 10

Paso 2

Definimos la rutina para simular el proceso de Hawkes mediante el método de descomposición eficiente.

simular_hawkes_exacto <- function(N, lambda_inicial, a, alpha, delta) {
  # 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)  
  
  # Condiciones iniciales según (1)
  T[1] <- 0
  lambda_antes[1] <- lambda_inicial
  lambda_despues[1] <- lambda_inicial
  
  for (i in 2:(N + 1)) {
    # a) 
    u0 <- runif(1)
    s0 <- - (1 / a) * log(u0)
    
    # b) 
    u1 <- runif(1)
    d <- 1 + (delta * log(u1)) / (lambda_despues[i - 1] - a)
    
    # c) 
    if (d > 0) {
      s1 <- -(1 / delta) * log(d)
      tau_i <- min(s0, s1)
    } else {
      tau_i <- s0
    }
    
    # d)
    T[i] <- T[i - 1] + tau_i
    
    # e) 
    lambda_antes[i] <- ((lambda_despues[i - 1] - a) * exp(-delta * tau_i)) + a
    
    # f) 
    lambda_despues[i] <- lambda_antes[i] + alpha
  }
 return(data.frame(
  i = 0:N,
  Tiempo = T
))

}

Paso 3

Realizar la simulación.

resultado <- simular_hawkes_exacto(N, lambda_inicial, a, alpha, delta)
resultado
##     i    Tiempo
## 1   0  0.000000
## 2   1  8.267558
## 3   2 14.084759
## 4   3 15.057884
## 5   4 40.002078
## 6   5 44.492300
## 7   6 55.263462
## 8   7 67.766349
## 9   8 71.210934
## 10  9 71.911811
## 11 10 91.692096

Paso 4

Graficamos el proceso de conteo.

eventos <- resultado$Tiempo  # vector de tiempos de los eventos
t_max <- max(eventos)
t <- seq(0, t_max, length.out = 500)  # grilla temporal
# Calcular el número acumulado de eventos hasta cada t
N_t <- sapply(t, function(x) sum(eventos <= x))

plot(t, N_t, type = "s",
     main = expression("Proceso de Hawkes"),
     xlab = expression(t),
     ylab = expression(N[t]),
     lwd = 2,
     col = "blue")

points(eventos, 0:(length(eventos) - 1), col = "red", pch = 19)

# Línea base gris
abline(h = 0, col = "gray")

Paso 5

Graficamos la función de intensidad

lambda_t <- sapply(t, function(x) {
  lambda_inicial + sum(alpha * exp(-delta * (x - eventos[eventos < x])))
})


plot(t, lambda_t, type = "l",
     main = expression("Intensidad del proceso de Hawkes"),
     xlab = expression(t),
     ylab = expression(lambda(t)),
     lwd = 2,
     col = "red")


abline(v = eventos, col = "gray80", lty = 3)

Algunos comentarios sobre el algoritmo

  • Los pasos del a., b. y c. se deduce directamente de la descomposición hecha anteriormente, y se aplica el proceso clásico de muestrear con la transformada inversa. Cada evento nuevo puede suceder en un tiempo \(s_0\) o \(s_1\) de manera independiente, y se simula cual sucede primero (el mínimo de los tiempos), de donde elegimos \(\tau_i\) para el paso d..

  • Los pasos e. y f. se deducen directamente de la definición de la función de intensidad \(\lambda(t)\), y estos corresponden a la intensidad justamente antes de que sucede el evento \(i\) y justamente después a que suceda el evento \(i\).

Referencias

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

[2] Rizoiu, M.-A., Lee, Y., Mishra, S., & Xie, L. (2018). A Tutorial on Hawkes Processes for Events in Social Media. arXiv preprint. https://arxiv.org/abs/1708.06401