# 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?

• I have had success fitting $\mu$ and $\alpha$ by maximizing the MLE the LBFGS implementation in scipy. The log-likelihood is not concave in $\beta$ though, so I simply iterated over a range of $\beta$ values and picked the one with the maximum likelihood. Note that $\alpha < \beta$ is required for stationarity of the process. Commented Jan 5, 2017 at 3:56
• curious, what is the correct form of the λ(t) function using the values of R(i) instead of resumming at each step?
– crow
Commented Sep 17, 2017 at 21:53

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 .

felix, using exp/log transforms seems to ensure positivity of the parameters. As for the small alpha thing, search the arxiv.org for a paper called "limit theorems for nearly unstable hawkes processes"

• 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.) Commented Jan 14, 2013 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. Commented Sep 23, 2012 at 15:26
• Could you give more details of what you did? It seems there is a problem with setting constraints and also that large beta is indistinguishable from zero alpha (they both look Poisson).
– Simd
Commented Aug 14, 2014 at 16:49

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)

• How do you ensure mu, alpha and beta are not set to negative values?
– Simd
Commented Aug 14, 2014 at 8:50
• You can set the lower and upper parameters in the optim call. Commented Aug 14, 2014 at 16:20
• Not for Nelder-Mead you can't can you which is the default? (See stat.ethz.ch/R-manual/R-devel/library/stats/html/optim.html ) . Also, I don't think there is any way to distinguish huge beta from zero alpha is there so a general optimization seems doomed.
– Simd
Commented Aug 14, 2014 at 16:44

Here is my solution to "What is the simplest practical method to implement?" using python, specifically numpy, scipy and tick.

One modification is that I set the exponential kernel such that alpha x beta x exp (-beta (t - ti)), to coincide with how tick defines exponential kernels: https://x-datainitiative.github.io/tick/modules/generated/tick.hawkes.HawkesExpKern.html#tick.hawkes.HawkesExpKern

I'm assuming the reader might not be familiar with python.

Import the relevant libraries:

import numpy as np
from scipy.optimize import minimize
from tick.hawkes import SimuHawkesExpKernels



Define the recursive function (R(i) in the question above) which returns an array the same size as the number of events:

def _recursive(timestamps, beta):
r_array = np.zeros(len(timestamps))
for i in range(1, len(timestamps)):
r_array[i] = np.exp(-beta * (timestamps[i] - timestamps[i - 1])) * (1 + r_array[i - 1])
return r_array


Define the log likelihood function specifying the various parameters:

def log_likelihood(timestamps, mu, alpha, beta, runtime):
r = _recursive(timestamps, beta)
return -runtime * mu + alpha * np.sum(np.exp(-beta * (runtime - timestamps)) - 1) + \
np.sum(np.log(mu + alpha * beta * r))


Simulate some Hawkes data using the tick library (or some other means):

m = 0.5
a = 0.2
b = 0.3
rt = 1000

simu = SimuHawkesExpKernels([[a]], b, [m], rt, seed=0)
simu.simulate()
t = simu.timestamps[0]


Scipy's minimize function is probably the most common python optimisation for scalar functions. It expects a function with only two sets of parameters; those you want to minimise and those which are fixed. There are various possible minimisation methods, I am using the default which is L-BFGS-B for bounded problems and is a quasi-newton Method.

Note that I should have started with a brute search first but the question asked for the simplest practical method. I could have also split into a minimisation over mu and alpha and then used simulated annealing over beta since the log-logarithm is convex over mu and alpha only.

Define a new function to be used by the minimize function and returns the negative log-likelihood:

def crit(params, *args):
mu, alpha, beta = params
timestamps, runtime = args
return -log_likelihood(timestamps, mu, alpha, beta, runtime)


Call the minimize function and set the bounds for m,a and b to be positive:

minimize(crit, [m + 0.1, a + 0.1, b+ 0.1], args=(t, rt), bounds=((1e-10, None), (1e-10, None), (1e-10, None),))


Which gives estimates of m,a,b: array([0.43835767, 0.25823306, 0.14769243])

• Although implementation is often mixed with substantive content in questions, we are supposed to be a site for providing information about statistics, machine learning, etc., not code. It can be good to provide code as well, but please elaborate your substantive answer in text for people who don't read this language well enough to recognize & extract the answer from the code. Commented Jul 2, 2020 at 18:07
• Hey @gung-ReinstateMonica thanks for the suggestion I have tried to put in some explanation please let me know what you think! Commented Jul 3, 2020 at 14:22