2  Univariate Hawkes

A univariate Hawkes process (Hawkes (1971)) is defined to be a self-exciting temporal point process where the conditional intensity function is given by

$$\lambda (t)=\mu (t)+{\mathrm{\Sigma }}_{i:{\tau }_i<t}\nu (t-{\tau }_i).$$

Here $\mu (t)$ is the background rate of the process and ${\mathrm{\Sigma }}_{i:{\tau }_i<t}\nu (t-{\tau }_i)$ is some historic temporal dependence (i.e., for times ${\tau }_i<t$, $ i = 1, …, T$). The classic homogeneous formulation uses an exponential decay kernal:

$$\begin{array}{}\text{(2.1)}& \lambda (t)=\mu +\alpha {\mathrm{\Sigma }}_{i:{\tau }_i<t}\text{exp}(-\beta \ast (t-{\tau }_i)).\end{array}$$

Here the parameter $\alpha$ is the increase in intensity immediately after the occurrence of an event, and $\beta >0$ controls the exponential decay of the intensity if no event has occurred. To avoid the conditional intensity heading off to infinity $\beta >\alpha$.

Plotted below is the conditional intensity ($\lambda (t),t\in [0,T]$) of a Hawkes process with $\mu =0.1$, $\alpha =1$, and $\beta =1.5$. Observed events (n = 13) are shown by the vertical dashes along the x-axis. The intensity increases immediately after an event occurs and decays exponentially over time if no event is observed for some period.

2.1 Simulating from a univariate Hawkes model

The function sim_hawkes() offers two simulation methods for simulating a realisation of a univariate Hawkes process. The default, with method = 1, uses the same method as algorithm 2 in Ogata (1981), simulating until a given time horizon (set by argument n).

sim <- sim_hawkes(mu = 0.5, alpha = 1, beta = 1.5, n = 40, plot = TRUE)

The second option, setting method = 2, uses an accept/reject framework and the argument n specifies the number of points to simulate (default n = 100).

sim <- sim_hawkes(mu = 0.5, alpha = 1, beta = 1.5, plot = TRUE, method = 2)

2.2 Fitting a univariate Hawkes model

require(stelfi)

To fit a univariate Hawkes model in stelfi use the function fit_hawkes() with the following required arguments

• times - a vector of numeric occurrence times, and
• parameters - a vector of named starting values for $\mu$ (mu), $\alpha$ (alpha), and $\beta$ (beta).

The function get_coefs() can then be called on the fitted model object to return the estimated parameter values.

2.2.1 A simulated example

Simulating a realisation of a univariate Hawkes process (see Section 2.1 for more details) with $\mu =1.3$, $\alpha =0.4$, and $\beta =1.5$ (over $t\in [0,500]$).

times <- sim_hawkes(mu = 1.3, alpha = 0.4, beta = 1.5, n = 500)

## starting values
sv <- c(mu = 1.3, alpha = 0.4, beta = 1.5)
## using stelfi
fit <- fit_hawkes(times = times, parameters = sv)
stelfi <- get_coefs(fit)
stelfi

       Estimate Std. Error
mu    1.3997320  0.1125668
alpha 0.3737118  0.1207772
beta  1.7730216  0.7495558

As a comparison, below emhawkes (Lee (2023)), and hawkesbow (Cheysson (2021)) are used to fit a univariate Hawkes process to the same simulated data.

## benchmark using emhawkes
require(emhawkes)
h <- new("hspec", mu = sv[1], alpha = sv[2], beta = sv[3])
## emhawkes requires the inter arrival times to fit the model
inter <- diff(times)
fit_em <- hfit(object = h, inter_arrival = inter)
em <- summary(fit_em)$estimate
em

        Estimate Std. error   t value      Pr(> t)
mu1    1.3976553  0.1091437 12.805644 1.524557e-37
alpha1 0.3683585  0.1380857  2.667608 7.639335e-03
beta1  1.7380101  0.8236957  2.110015 3.485708e-02

## bench mark using hawkesbow
require(hawkesbow)
fit_bow <- mle(events = times, kern = "Exponential", end = max(times))
## use the Hessian to obtain the standard errors from hawkesbow
bow <- cbind(Estimate = fit_bow$par,
         "Std. error" = -fit_bow$model$ddloglik(times, max(times)) |> solve() |> diag() |> sqrt())
bow

      Estimate Std. error
[1,] 1.3998368 0.11238176
[2,] 0.2107781 0.05814911
[3,] 1.7724499 0.74877227

The table below gives the estimated parameter values from each of stelfi, emhawkes, and hawkesbow along with the standard errors in brackets. Note that hawkesbow estimates $\frac{\alpha }{\beta }$ rather that $\beta$ directly, and that the standard errors are computed here using the returned Hessian matrix $H$ (i.e., $\sqrt{\text{diag}-(H^{-1})}$).

	$\mu$	$\alpha$	$\beta$	$\frac{\alpha }{\beta }$
TRUTH	1.300	0.400	1.500	0.267
stelfi	1.4 ( 0.113 )	0.374 ( 0.121 )	1.773 ( 0.75 )	-
emhawkes	1.398 ( 0.109 )	0.368 ( 0.138 )	1.738 ( 0.824 )	-
hawkesbow	1.4 ( 0.112 )	-	1.772 ( 0.749 )	0.211 ( 0.058 )

2.2.2 An applied example

A NIWA scientist found a working USB in the scat of a leopard seal, they then tweeted about it in the hopes of finding its owner¹. The dates and times of these tweets and retweets are available in stelfi as retweets_niwa.

data(retweets_niwa)

The dates/times need to be numeric and sorted in ascending order (starting a time $t=0$). Note too that there can be no simultaneous events.

## numeric time stamps
times <- unique(sort(as.numeric(difftime(retweets_niwa ,min(retweets_niwa),units = "mins"))))

The histogram below shows the observed counts (4772) of unique retweet times from the original tweet ($t=0$) to the time (just over two days later) that the owner of the USB came forward, t = 3143 mins.

To fit the model chose some starting values for the parameters and supply these along with the numeric times to the function fit_hawkes().

sv <- c(mu = 9, alpha = 3, beta = 10)
fit <- fit_hawkes(times = times, parameters = sv) 

## print out estimated parameters
pars <- get_coefs(fit)
pars

        Estimate  Std. Error
mu    0.06328099 0.017783908
alpha 0.07596531 0.007777899
beta  0.07911346 0.008109789

From the estimated coefficients above

• the expected (estimated) background rate of sightings (i.e., independent sightings) is $\hat{\mu }T=$ 0.063 $\times$ 3143.02 $=$ 198.89, which indicates that ~199 retweets (from the 4772) were principal retweets and that the remaining were due to self-excitement;
• the expected number of retweets “triggered” by any one retweet is estimated as $\frac{\hat{\alpha }}{\hat{\beta }}$ = 0.96 (note the maximum this can possibly be is 1);
• the expected number of descendants per retweet is estimated as $\frac{\hat{\beta }}{\hat{\beta }-\hat{\alpha }}$ = 25.13;
• the rate of decay for the self-excitement is estimated as $\frac{1}{\hat{\beta }}$ = 12.64, indicating that after ~13 minutes a retweet is likely unrelated² to the previous retweet.

The show_hawkes() function can be called on the fitted model object to plot the estimated conditional intensity and the data, top and bottom panels below respectively.

show_hawkes(fit)

2.3 Goodness-of-fit for a univariate Hawkes process

The compensator, $\mathrm{\Lambda }(\tilde{t})$, of any inhomogeneous Poisson process gives the expected number of events in some defined interval $[0,\tilde{t}]$ is given by

$$\mathrm{\Lambda }(\tilde{t})={\int }_0^{\tilde{t}}\lambda (t)dt.$$

The random change theorem (Daley, Vere-Jones, et al. (2003)) states that if a set of events $[{\tau }_1,...,{\tau }_n]$ is a realisation from a inhomogeneous Poisson process then $[\mathrm{\Lambda }({\tau }_1),...,\mathrm{\Lambda }({\tau }_n)]$ is a realisation of a homogeneous Poisson process with unit rate. Letting $\delta {\mathrm{\Lambda }}_i=\mathrm{\Lambda }({\tau }_i)-\mathrm{\Lambda }({\tau }_{i-1})$ for $i=2,...,N$ and $\delta {\mathrm{\Lambda }}_1=\mathrm{\Lambda }({\tau }_1)$, under the theorem above $\delta {\mathrm{\Lambda }}_i\sim \text{Exp}(1)$. Using this result a typical goodness-of-fit test is a Kolmogorov-Smirnov (KS) test (see Daley, Vere-Jones, et al. (2003) for more details) where the KS statistic is a measure of the distance between the empirical distribution of all $\delta {\mathrm{\Lambda }}_i$ and the CDF of a $\text{Exp}(1)$ distribution.

The compensator differences can be extracted from a fitted model using the compensator_differences() function, and a KS test can be carried out manually:

times <- sim_hawkes(mu = 1.3, alpha = 0.4, beta = 1.5, n = 500)
sv <- c(mu = 2, alpha = 1, beta = 5)
fit <- fit_hawkes(times = times, parameters = sv)
compensator <- compensator_differences(fit)
stats::ks.test(compensator, "pexp")

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  compensator
D = 0.013855, p-value = 0.9958
alternative hypothesis: two-sided

This gives no evidence against the compensator values coming from a $\text{Exp}(1)$ distribution.

Another goodness-of-fit test is the Box-Ljung (or Ljung–Box) test, which tests for autocorrelation between the consecutive compensator values (i.e., independence/stationarity).

stats::Box.test(compensator, type = "Ljung")

    Box-Ljung test

data:  compensator
X-squared = 0.66599, df = 1, p-value = 0.4145

This gives no evidence against the consecutive compensator values being independently distributed.

Alternatively, both tests and some diagnostic plots are returned by calling the show_hawkes_GOF() function. The four panels an be interpreted as follows

• top left, plots the compensator values against the observed times, which under a well fitting model should align;
• top right, a transformed QQ plot, the observed quantities should align with the theoretical quantiles under a well fitting model;
• bottom left, the compensator differences, which under the model are assumed to be $\text{Exp}(1)$ distributed;
• bottom right, consecutive compensator differences, which should show no obvious pattern (no autocorrelation evident) under a well fitting model.

show_hawkes_GOF(fit)

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  interarrivals
D = 0.013855, p-value = 0.9958
alternative hypothesis: two-sided


    Box-Ljung test

data:  interarrivals
X-squared = 0.66599, df = 1, p-value = 0.4145

2.4 Fitting an inhomogenous Hawkes model

A univariate inhomogenous Hawkes process has $\mu =\mu (t)$ in Equation 2.1 (i.e., the baseline is time varying, rather than being constant).

Below we simulate data from a Hawkes process with $\mu (t)=[A+B\ast \text{sin}(\frac{2\pi t}{365.25})]$ (with $A=1$ and $B=0.5$) to represent a yearly cycle and self-exciting parameters $\alpha =1$ and $\beta =2$ using hawkesbow (Cheysson (2021)).

mut <- function(t) {
 1 + 0.5*sin((2*pi*t)/365.25)
}

set.seed(1)
times <- hawkesbow::hawkes(1000, fun = mut,
 M = 1.5, repr = 0.5, family = "exp", rate = 2)$p

To fit this model in stelfi the function $\mu (t)$ (background) and its integral ${\int }_0^t\mu (y)dy$ (background_integral) are supplied by the user. To ensure $\mu (t)>0$ below $\mu (t)$ is written as a function of $x$ and $y$, where $A=e^x$ and $B=\text{logit}(y)e^x$.

background <- function(params,t){
        A = exp(params[[1]])
        B = stats::plogis(params[[2]]) * A
        return(A + B*sin((2*pi*t)/365.25))
}

background_integral <- function(params,x){
        A = exp(params[[1]])
        B = stats::plogis(params[[2]]) * A
        return((A*x)-B*cos((2*pi*x)/365.25))
}

These functions are then passed to the function fit_hawkes_cbf() alongside the observed times (times) and a list of starting values for both $\mu (t)$ and the self-exciting components via the argument background_parameters and parameters respectively.

sv = list(alpha = 0.5, beta = 1.5)
background_sv = list(1,1)
fit <- fit_hawkes_cbf(times = times, parameters = sv,
background = background, background_integral = background_integral, background_parameters = background_sv)

The estimated values of the the transformed $A$ and $B$ parameters are returned via fit$background_parameters, the code below transforms them back to the original scale.

exp(fit$background_parameters[1])

[1] 1.135523

plogis(fit$background_parameters[2]) * exp(fit$background_parameters[1])

[1] 0.7755937

The estimated values of the self-excitement parameters $\alpha$ and $\beta$ are returned by the get_coefs() function as usual. The table below compares the estimated values to the true ones

	$A$	$B$	$\alpha$	$\beta$
TRUTH	1.000	0.500	1.000	2.000
stelfi	1.136	0.776	0.985	2.198

The function show_hawkes() can be used to show the fitted model.

show_hawkes(fit)

The function show_hawkes_GOF() will show/print the model diagnostics discussed in Section 2.3.

show_hawkes_GOF(fit)

    Asymptotic one-sample Kolmogorov-Smirnov test

data:  interarrivals
D = 0.018277, p-value = 0.4985
alternative hypothesis: two-sided


    Box-Ljung test

data:  interarrivals
X-squared = 28.608, df = 1, p-value = 8.86e-08

Cheysson, Felix. 2021. Hawkesbow: Estimation of Hawkes Processes from Binned Observations. https://CRAN.R-project.org/package=hawkesbow.

Daley, Daryl J, David Vere-Jones, et al. 2003. An Introduction to the Theory of Point Processes: Volume i: Elementary Theory and Methods. Springer.

Hawkes, AG. 1971. “Spectra of Some Self-Exciting and Mutually Exciting Point Processes.” Biometrika.

Lee, Kyungsub. 2023. Emhawkes: Exponential Multivariate Hawkes Model. https://CRAN.R-project.org/package=emhawkes.

Ogata, Y. 1981. “On Lewis’ Simulation Method for Point Processes.” IEEE Transactions on Information Theory. https://doi.org/10.1109/TIT.1981.1056305.

---

1. This was at the time it was called Twitter (not X) so ‘tweet’ is the correct term!

2. Note that unrelated is not quite the correct term here, rather that the event is likely a baseline event rather than due to self-excitement