The thinning algorithm

  1. Dado un proceso de Hawkes de la forma: \[ \lambda(t \mid \mathcal{H}_t) = \lambda_0(t) + \sum_{i: \, t > T_i} \phi(t - T_i) \tag{1} \]

  2. Fijar el tiempo actual \(T = 0\) y el contador de eventos \(i = 1\).

  3. Mientras \(i \le N\)

    1. Calcular la cota superior de la intensidad de Poisson: (Usando Ec(1)) \[ \lambda^{*} = \lambda(T) \]

    2. Muestrear tiempo entre llegadas: Tome un \(u \sim U(0,1)\) y sea \(\tau=-\frac{ln(u)}{\lambda^{*}}\)

    3. Actualizar el tiempo presente \(T=T+\tau\)

    4. Tome \(s \sim U(0,1)\)

    5. Si \(s\leq \frac{\lambda(T)}{\lambda^*}\) acepte la muestra actual y sea \(T_i=T\) y \(i=i+1\). Otro caso rechaze la muestra y vuelva al paso a.

Paso 1

Por simplicidad se usa el kernel exponencial definido de la siguiente forma: \[ \phi(x)=\alpha e^{-\delta x} \] Donde \(\delta\) y \(\alpha\) son hiperparametros, en este caso se toman como \(\delta=0.6\) y \(\alpha=1.2\). Se puede observar que es una función monotona decreciente, en la siguiente figura:

# Definición del kernel
phi <- function(x, alpha, delta) {
  alpha * exp(-delta * x)
}

# Parámetros
alpha <- 0.6
delta <- 1.2

# Gráfico
curve(phi(x, alpha, delta),
      from = 0, to = 10,
      xlab = expression(x),
      ylab = expression(phi(x)),
      main = "Kernel Exponencial",
      lwd = 2, col = "blue")

abline(h = 0, col = "gray")

Por ultimo para tener el procesos totalmente definido se necesita el \(\lambda_0(t)\), por practicidad tambien se tomara como una constante.\(\lambda_0(t)=0.2\)

Con esto nuestro proceso de Hawkes esta completamente definido por su función de intensidad de la forma:

\[ \lambda(t \mid \mathcal{H}_t) = 0.2 + \sum_{i: \, t > T_i} 0.6 e^{-1.2 (t - T_i)} \tag{2} \]

Paso 2

Se crean objetos que guarden el tiempo \(T=0\) y el contador de evento \(i=1\)

T=0
i=1

Paso 3

\(N\) es el numero de eventos que el agoritmo va a simular. Por facilidad de va a fijar un \(N=10\), es decir solo se van a simular la ocurrencia de 10 eventos.

N=10

a.

La cota superior se toma como \(\lambda^*=\lambda(T)\), es decir se toma como cota, el valor de la intensidad de cuando ocurre el evento justamente anterior, en este caso se toma la intensidad base, es decir \(\lambda^*=\lambda_0(t)=0.2\)

lambda_est = 0.2

b.

Se muestrea un valor de una distribución uniforme y se calcula el tiempo entre llegadas.

u=runif(1,0,1)
tau=-log(u)/lambda_est

c.

Se actualiza el tiempo \(T\), con el \(\tau\) encontrado en el paso anterior.

T=T+tau

d.

Se halla s muestreando de una distribución uniforme.

s=runif(1,0,1)

e.

se mira la desigualdad, si es verdadera se acepta el valor. Como estamos en el primero paso, la función de intesidad concuerda con la cota superioror, es decir \(\lambda^*=\lambda(T)\) ya que la suma \(\sum_{i:t>T_i} 0.6 e^{-1.2(t-T_i)}=0\), es vacia, por lo tanto es 0. (Preguntar, si seria \(\lambda(T)=\lambda_0+0.6\))

s <= 1
## [1] TRUE
T=T
i=i+1
print(T)
## [1] 2.158238
print(i)
## [1] 2

Entonces el primer evento del proceso, ocurre en el nuevo tiempo T. y ya empieza el contador para el nuevo evento con i=2

Simulación Completa

Para hacer la simulación de una manera eficiente se crea un función que realiza el algoritmo de manera automatica y guarda los tiempos en el que sucede cada evento.

inten <- function(t, eventos, alpha, delta, lambda_0) {
  if (length(eventos) == 0) {
    return(lambda_0)
  } else {
    return(lambda_0 + sum(alpha * exp(-delta * (t - eventos[eventos < t]))))
  }
}



SimHawkes <- function(N, alpha, delta, lambda_0) {
  eventos <- numeric(0)  
  T <- 0                 
  
  while (length(eventos) < N) {
    # (a) 
    lambda_ast <- inten(T, eventos, alpha, delta, lambda_0)
    
    # (b) 
    u <- runif(1)
    tau <- -log(u) / lambda_ast
    
    # (c)
    T <- T + tau
    
    # (d) 
    lambda_T <- inten(T, eventos, alpha, delta, lambda_0)
    
    # (e) 
    s <- runif(1)
    if (s <= lambda_T / lambda_ast) {
      
      eventos <- c(eventos, T)
     
    } else {
      
      next
    }
  }
  
  return(eventos)
}

Ya con la función es correrla para nuestro caso particular

set.seed(2)
eventos=SimHawkes(N = 10, alpha = 0.6, delta = 1.2, lambda_0 = 0.2)
print(eventos)
##  [1]  8.44018 11.22168 11.48285 17.99094 20.94955 22.20972 25.00148 28.58887
##  [9] 30.56845 31.26247

Adicionalmente se pueden graficar, los procesos y su función de intensidad.

t_max <- max(eventos)
t <- seq(0, t_max, length.out = 500) #crear grilla del 0 al maximo con 500 puntos

N_t <- sapply(t, function(x) sum(eventos <= x)) #Calcula la suma acumulada, para cada elemnto del vector t

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

# Agregar los puntos de evento en rojo
points(eventos, 1:length(eventos), col = "red", pch = 19)

# Línea horizontal gris de referencia
abline(h = 0, col = "gray")

t_max <- max(eventos)+5
t <- seq(0, t_max+5, length.out = 500)

lambda_t <- sapply(t, function(x) inten(x, eventos,alpha = 0.6, delta = 1.2, lambda_0 = 0.2))


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

Referencias

Rizoiu, Marian-Andrei; Lee, Young; Mishra, Swapnil; Xie, Lexing (2017). A Tutorial on Hawkes Processes for Events in Social Media. Frontiers of Multimedia Research, pp. 191-218. DOI: 10.1145/3122865.3122874.