Solved – Finding the MLE for a univariate exponential Hawkes process

likelihoodmaximum likelihoodstochastic-processes

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?

Best Answer

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"

Related Solutions

Solved – MLE of a multivariate Hawkes process

There is a small mistake in the derivation. In line 5 (in the inserted figure) one needs $T = t_{1,F} = t_{2,G}$ for the identity to be correct, and this is generally not the case. The terms in the final sums should be $e^{-\beta_{i,1}(T - t_{1,f})} - 1$ and $e^{-\beta_{i,2}(T - t_{2,g})} - 1$, respectively. Otherwise the derivation looks correct.

A slightly simpler derivation can take line 3 as a starting point. Then interchange the sums and integration with the resulting inner integral being from $t_{j,k}$ to $T$.

It might be worth noting that for the Hawkes process considered here, it is possible to compute $\lambda_i^*(t_{i,j})$ recursively, which implies that the computational complexity of the log-likelihood can be made linear in the number of jumps (instead of quadratic as the double sum over the jumps suggests).

I doubt that there are inconsistent versions of the likelihood in the literature, but there may, of course, be mistakes in some of the references. Another (likely) possibility is that the notation or the assumptions differ, or that the representations are, indeed, equivalent, but written in different ways. One possibility is that the baseline intensity $\lambda_i$ is omitted, so that the $\lambda_i T$ term disappears.

Solved – How to compute partial log-likelihood function in Cox proportional hazards model

This is technically a programming question with an easy programming answer. If you simply want the partial likelihood, why not fool R into giving it to you? Simply initialize beta and allow no iterations, then extract the loglik value from the coxph object. (see ?coxph.object).

For example:

## artificial data
library(survival)
n <- 1000
t <- rexp(100)
c <- rbinom(100, 1, .2) ## censoring indicator (independent process)
x <- rbinom(100, 1, exp(-t)) ## some arbitrary relationship btn x and t
betamax <- coxph(Surv(t, c) ~ x)
beta1 <- coxph(Surv(t, c) ~ x, init = c(1), control=list('iter.max'=0))

With example output:

> betamax$loglik
[1] -68.62548 -65.99652
> beta1$loglik
[1] -66.10908 -66.10908

You can even define a wrapper:

loglik <- function(beta, formula) {
  formula, init=beta, control=list('iter.max'=0))$loglik[2]
}

betas <- seq(0, 2, by=0.01)
logliks <- sapply(betas, loglik, Surv(t, c) ~ x)
plot(betas, logliks)
abline(v=betamax$coefficients)

Best Answer

Related Solutions

Solved – MLE of a multivariate Hawkes process

Solved – How to compute partial log-likelihood function in Cox proportional hazards model

Related Question