try package RHawkes

RHawkes.R

library(RHawkes)

simulate_branching_death <- function(mu,
                                     alpha,
                                     beta,
                                     Tmax) {
  
  N0 <- rpois(1, Tmax * mu)
  if(N0==0L){
    return(
      data.frame(
        id = NA,
        parent = 0L,
        birth_time = NA,
        death_time = NA,
        generation = NA
      )
    )
  }
  
  t0 <- sort(runif(N0,0,Tmax))
  individuals <- data.frame(
    id = seq_len(N0),
    parent = 0L,
    birth_time = t0,
    death_time = t0 + rexp(N0, rate = beta),
    generation = 1L
  )
  
  next_id <- N0 + 1L
  i <- 1L
  
  while (i <= nrow(individuals)) {
    
    birth <- individuals$birth_time[i]
    death <- min(individuals$death_time[i], Tmax)
    
    if (birth < death) {
      
      t <- birth
      
      repeat {
        
        t <- t + rexp(1, rate = alpha)
        
        if (t > death || t > Tmax) {
          break
        }
        
        individuals <- rbind(
          individuals,
          data.frame(
            id = next_id,
            parent = individuals$id[i],
            birth_time = t,
            death_time = t + rexp(1, rate = beta),
            generation =
              individuals$generation[i] + 1L
          )
        )
        
        next_id <- next_id + 1L
      }
    }
    
    i <- i + 1L
  }
  return(individuals[order(individuals$birth_time),])
}

Tmax <- 10
set.seed(123456)
tree <- simulate_branching_death(
  mu = 1,
  alpha = 1,
  beta = 1,
  Tmax = Tmax
)

system.time( est1 <- EM1partial(tree$birth_time, Tmax, c(2,1,0.5,1)) )

plot_cumcount <- function(ti, Tmax){
  plot(c(0,ti, Tmax), c(0,seq_along(ti), length(ti)), type = "s",
       xlab="time", ylab="cumulative count")
}


cumint_hawkes <- function(xv, mu, a, par_wei, ti){
  sapply(xv, function(t){mu*t+sum(a*pweibull(t-ti[ti<t], shape = par_wei[1], scale = par_wei[2]))})
}

xv <- seq(0, Tmax, by=0.1)
Lvhat <- cumint_hawkes(xv,est1$pars[1],est1$pars[4],est1$pars[2:3],tree$birth_time)
Lv <- cumint_hawkes(xv,1,1,c(1,1),tree$birth_time)
png("step.png")
plot_cumcount(tree$birth_time, Tmax)
lines(xv,Lv,lty=2, col="royalblue")
lines(xv,Lvhat,lty=3, col="orangered")
legend("topleft", c("true","MLE"), col=c("royalblue","orangered"), lty=2:3)
dev.off()