Finding the MLE for a univariate exponential Hawkes process

The univariate exponential Hawkes process is a self-exciting point process with an event arrival rate of:

$\lambda(t) = \mu + \sum\limits_{t_i<t}{\alpha e^{-\beta(t-t_i)}}$

where $t_1,..t_n$ are the event arrival times.

The log likelihood function is

$- t_n \mu + \frac{\alpha}{\beta} \sum{( e^{-\beta(t_n-t_i)}-1 )} + \sum\limits_{i<j}{\ln(\mu+\alpha e^{-\beta(t_j-t_i)})}$

which can be calculated recursively:

$- t_n \mu + \frac{\alpha}{\beta} \sum{( e^{-\beta(t_n-t_i)}-1 )} + \sum{\ln(\mu+\alpha R(i))}$

$R(i) = e^{-\beta(t_i-t_{i-1})} (1+R(i-1))$

$R(1) = 0$

What numerical methods can I use to find the MLE? What is the simplest practical method to implement?

-

The Nelder-Mead simplex algorithm seems to work well.. It is implemented in Java by the Apache Commons Math library at https://commons.apache.org/math/ . I've also written a paper about the Hawkes processes at Point Process Models for Multivariate High-Frequency Irregularly Spaced Data . On an unrelated note, do you know if the elegant recursions available for the exponential kernel are possible with other kernels, or is this a property exclusively pertaining to the exponential kernel?

-
Welcome to our site! –  whuber Jan 14 '13 at 18:55
Welcome to the site, @StephenCrowley. If you have your own question, please do not post it as ( / as part of) an answer. Click on the gray "ASK QUESTION" button at the top of the page & ask it there. If you have a question for clarification from the OP, you should ask it in a comment to the question post above. (Although frustratingly, you cannot do that until you reach 50 rep.) –  gung Jan 14 '13 at 18:58

I solved this problem using the nlopt library. I found a number of the methods converged quite quickly.

-
I assume you're familiar with T. Ozaki (1979), Maximum likelihood estimation of Hawkes' self-exciting point processes, Ann. Inst. Statist. Math., vol. 31, no. 1, 145-155. –  cardinal Sep 23 '12 at 15:26

You could also do a simple maximization. In R:

neg.loglik <- function(params, data, opt=TRUE) {
mu <- params[1]
alpha <- params[2]
beta <- params[3]
t <- sort(data)
r <- rep(0,length(t))
for(i in 2:length(t)) {
r[i] <- exp(-beta*(t[i]-t[i-1]))*(1+r[i-1])
}
loglik <- -tail(t,1)*mu
loglik <- loglik+alpha/beta*sum(exp(-beta*(tail(t,1)-t))-1)
loglik <- loglik+sum(log(mu+alpha*r))
if(!opt) {
return(list(negloglik=-loglik, mu=mu, alpha=alpha, beta=beta, t=t,
r=r))
}
else {
return(-loglik)
}
}

# insert your values for (mu, alpha, beta) in par
# insert your times for data
opt <- optim(par=c(1,2,3), fn=neg.loglik, data=data)

-