A simulation of maximum likelihood estimation of Hawkes process

sim_Hawkes_MLE.R

library(hawkes)
mu <- 0.5
a  <- 0.7
b  <- 1.0
Tau <- 100
set.seed(12345)
dat <- simulateHawkes(
  lambda0 = mu,
  alpha = a*b,
  beta = b,
  horizon = Tau
)

ll <- function(par, ti, Tau){
  mu <- exp(par[1])
  b <- exp(par[2])
  a <- exp(par[3])
  sum(log(mu+sapply(ti, function(t)sum(a*b*exp(-b*(t-ti[ti<t])))))) - 
    (mu*Tau+sum(a*(1-exp(-b*(Tau-ti)))))
}


opt <- optim(rep(0,3), ll,
             ti = dat[[1]], Tau = Tau,
             control = list(fnscale=-1), method = "BFGS")
print(opt)

muhat <- exp(opt$par[1])
bhat <- exp(opt$par[2])
ahat <- exp(opt$par[3])

ti <- dat[[1]]

# lv <- mu+sapply(xv, function(t)sum(a*b*exp(-b*(t-ti[ti<t]))))
xv <- seq(0,Tau,by=1)
cumint_hawkes <- function(xv, mu, a, b, ti){
  sapply(xv, function(t){mu*t+sum(a*(1-exp(-b*(t-ti[ti<t]))))})  
}

Lv <- cumint_hawkes(xv,mu,a,b,ti)
Lvhat <- cumint_hawkes(xv,muhat,ahat,bhat,ti)

png("fit.png")
plot(c(0, ti), c(0L, seq_along(ti)), type="s", xlab = "time", ylab="couunt")
lines(xv, Lv, col="royalblue", lty=2)
lines(xv, Lvhat, col="orangered", lty=3)
legend("topleft",legend = c("observed","true","MLE"), 
       lty=1:3, col=c("black","royalblue","orangered"))
dev.off()