Simulating Ogata's Thinning Algorithm in R

Question

I am trying to implement Ogata's Thinning Algorithm exactly as given in Algorithm 3 in https://www.math.fsu.edu/~ychen/research/Thinning%20algorithm.pdf with the parameters they specify to generate Figure 2.

This is the code that replicates it verbatim.

# Simulation of a Univariate Hawkes Poisson with
# Exponential Kernel γ(u) = αe^(−βu), on [0, T].
# Based on: https://www.math.fsu.edu/~ychen/research/Thinning%20algorithm.pdf

library(tidyverse)

# Initialize parameters that remains the same for all realizations
mu <- 1.2
alpha <- 0.6
beta <- 0.8
T <- 10

# Initialize other variables that change through realizations
bigTau <- vector('numeric')
bigTau <- c(bigTau,0)
s <- 0
n <- 0
lambda_vec_accepted <- c(mu)

# Begin loop
while (s < T) {

  # -------------------------------
  # Compute lambda_bar
  # -------------------------------
  sum_over_big_Tau <- 0
  for (i in c(1:length(bigTau))) {
    sum_over_big_Tau <- sum_over_big_Tau + alpha*exp(-beta*(s-bigTau[i]))
  }
  lambda_bar <- mu + sum_over_big_Tau

  u <- runif(1)
  # so that w ∼ exponential(λ_bar)
  w <- -log(u)/lambda_bar
  # so that s is the next candidate point
  s <- s+w
  # Generate D ∼ uniform(0,1)
  D <- runif(1)

  # -------------------------------
  # Compute lambda_s
  # -------------------------------
  sum_over_big_Tau <- 0
  for (i in c(1:length(bigTau))) {
    sum_over_big_Tau <- sum_over_big_Tau + alpha*exp(-beta*(s-bigTau[i]))
  }
  lambda_s <- mu + sum_over_big_Tau

  # -------------------------------
  # Accepting with prob. λ_s/λ_bar
  # -------------------------------
  if (D*lambda_bar <= lambda_s ) {

    lambda_vec_accepted <- c(lambda_vec_accepted,lambda_s)

    # updating the number of points accepted
    n <- n+1
    # naming it t_n
    t_n <- s
    # adding t_n to the ordered set bigTau
    bigTau <- c(bigTau,t_n)
  }

}

df<-data.frame(x=bigTau,y=lambda_vec_accepted)
ggplot(df) + geom_line(aes(x=bigTau,y=lambda_vec_accepted))

However, the plot I managed to get (running several times) for lambda vs time is something like this and nowhere near what they got in Figure 2 (exponentially decreasing).

Valid XHTML

I am not sure what is the mistake I am doing. It will be great if anyone can help. This is needed for my research. I am from biology and so please excuse if I am doing something silly. Thanks.

Your code seems correct but there is a "resolution" problem. Indeed you store the value of the intensity in `bigTau` only when you keep an event. You should rather compute it and store it at a fixed time step `dt`. If you are familiar with Python, you can also try the [tick](https://x-datainitiative.github.io/tick/) library to simulate Hawkes processes and obtain a similar graph. — user2015762, Aug 22 '18 at 10:18

Simulating Ogata's Thinning Algorithm in R

0 Answers0