Encontrar el MLE para un proceso de Hawkes exponencial univariante

Question

El proceso exponencial univariante de Hawkes es un proceso puntual autoexcitado con una tasa de llegada de eventos de:

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

donde $ t_1,..t_n $ son los tiempos de llegada del evento.

La función de probabilidad logarítmica es

$ - 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)})} $

que se puede calcular de forma recursiva:

$ - 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 $

¿Qué métodos numéricos puedo utilizar para encontrar la MLE? ¿Cuál es el método práctico más sencillo de aplicar?

Answer 1

El algoritmo simplex de Nelder-Mead parece funcionar bien.. Está implementado en Java por la biblioteca Apache Commons Math en https://commons.apache.org/math/ . También he escrito un artículo sobre los procesos de Hawkes en Modelos de procesos puntuales para datos multivariantes de alta frecuencia y espaciados irregularmente .

felix, el uso de transformaciones exp/log parece garantizar la positividad de los parámetros. En cuanto a lo de la alfa pequeña, busca en arxiv.org un artículo llamado "limit theorems for nearly unstable hawkes processes"

Answer 2

He resuelto este problema utilizando el nlopt biblioteca. He comprobado que varios de los métodos convergen con bastante rapidez.

Answer 3

También puedes hacer una simple maximización. En 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)

Answer 4

Aquí está mi solución a "¿Cuál es el método práctico más simple de implementar?" usando python, específicamente numpy, scipy y tick.

Una modificación es que establezco el núcleo exponencial tal que alfa x beta x exp (-beta (t - ti)), para que coincida con la forma en que tick define los núcleos exponenciales: https://x-datainitiative.github.io/tick/modules/generated/tick.hawkes.HawkesExpKern.html#tick.hawkes.HawkesExpKern

Asumo que el lector puede no estar familiarizado con python.

Importe las bibliotecas correspondientes:

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

Defina la función recursiva (R(i) en la pregunta anterior) que devuelve una matriz del mismo tamaño que el número de eventos:

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

Defina la función de probabilidad logarítmica especificando los distintos parámetros:

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))

Simule algunos datos de Hawkes utilizando la biblioteca de ticks (o algún otro medio):

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]

La función minimizar de Scipy es probablemente la optimización más común de python para las funciones escalares. Espera una función con sólo dos conjuntos de parámetros; los que quieres minimizar y los que son fijos. Hay varios métodos de minimización posibles, yo estoy usando el predeterminado que es L-BFGS-B para problemas acotados y es un método cuasi-newton.

Ten en cuenta que debería haber empezado con una búsqueda bruta primero, pero la pregunta pedía el método práctico más sencillo. También podría haber dividido en una minimización sobre mu y alfa y luego utilizar el recocido simulado sobre beta ya que el logaritmo es convexo sobre mu y alfa solamente.

Define una nueva función que será utilizada por la función minimizar y devuelve la log-verosimilitud negativa:

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

Llama a la función minimizar y establece que los límites de m,a y b sean positivos:

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

Lo que da estimaciones de m,a,b: array([0.43835767, 0.25823306, 0.14769243])