In cmjt/stelfi: Hawkes and Log-Gaussian Cox Point Processes Using Template Model Builder

Hawkes process {#hawkes}

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

$$\lambda(t) = \mu(t) + \Sigma_{i:\tau_i<t}\nu(t-\tau_i)$$ where where $\mu(t)$ is the background rate of the process and $\Sigma_{i:\tau_i<t}\nu(t-\tau_i)$ is some historic temporal dependence. First introduced by @hawkes, the classic homogeneous formulation is:

$$\lambda(t) = \mu + \alpha \Sigma_{i:\tau_i<t}\text{exp}(-\beta * (t-\tau_i)) $$

The fit_hawkes() function

library(stelfi)
args(fit_hawkes)

Fitting a Hawkes model

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.

data(retweets_niwa)
head(retweets_niwa)

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

# NIWA seal first tweet according to Twitter
start <- lubridate::ymd_hms("2019-02-05 02:27:07")
## NIWA seal found owner tweet according to Twitter
end <- lubridate::ymd_hms("2019-02-07 06:50:08")
## hist
hist(retweets_niwa, breaks = "hours", axes = FALSE, 
     xlab = "", ylab = "",main = "", col = "grey",border = "grey",freq = TRUE)
axis(2,at = c(0,250))
mtext(2, line = 0, text = "Number of retweets per hour",cex = 0.8)
mtext(1,line = 1, at = c(start,end),text = c("NIWA \n tweeted","USB owner \n found"),cex = 0.9)

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

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

show_hawkes(fit)

show_hawkes_GOF(fit)

Fitting an ETAS-type marked model

Here we fit a univariate marked Hawkes process where the conditional intensity function is given by

$$\lambda(t; m(t)) = \mu + \alpha \Sigma_{i:\tau_i<t}m(\tau_i)\text{exp}(-\beta * (t-\tau_i)) $$ where $\mu$ is the background rate of the process and $m(t)$ is the temporal mark. Each event $i$ has an associated mark $\tau_i$ that multiples the self-exciting component of $\lambda$.

In this example, the events are earthquakes and the marks are the Richter magnitude of each earthquake.

data("nz_earthquakes")
head(nz_earthquakes)
nz_earthquakes <- nz_earthquakes[order(nz_earthquakes$origintime),]
nz_earthquakes <- nz_earthquakes[!duplicated(nz_earthquakes$origintime),]
times <- nz_earthquakes$origintime
times <- as.numeric(difftime(times , min(times), units = "mins"))
marks <- nz_earthquakes$magnitude
params <- c(mu = 3, alpha = 0.05, beta = 1)
fit <- fit_hawkes(times = times, parameters = params, marks = marks) 
## print out estimated parameters
pars <- get_coefs(fit)
pars

show_hawkes(fit)

show_hawkes_GOF(fit)

The fit_mhawkes() function

args(fit_mhawkes)

A multivariate Hawkes process allows for between- and within-stream self-excitement. In stelfi the conditional intensity for the $j^{th}$ ($j = 1, ..., N$) stream is given by

$$\lambda(t)^{j*} = \mu_j + \Sigma_{k = 1}^N\Sigma_{i:\tau_i<t} \alpha_{jk} e^{(-\beta_j * (t-\tau_i))},$$ where $j, k \in (1, ..., N)$. Here, $\alpha_{jk}$ is the excitement caused by the $k^{th}$ stream on the $j^{th}$. Therefore, $\boldsymbol{\alpha}$ is an $N \times N$ matrix where the diagonals represent the within-stream excitement and the off-diagonals represent the excitement between streams.

data(multi_hawkes)

fit <- stelfi::fit_mhawkes(times = multi_hawkes$times, stream = multi_hawkes$stream,
                           parameters = list(mu =  c(0.2,0.2),
                                        alpha =  matrix(c(0.5,0.1,0.1,0.5),byrow = TRUE,nrow = 2),
                                        beta = c(0.7,0.7)))
get_coefs(fit)

The fit_hawkes_cbf() function

args(fit_hawkes_cbf)

Fitting an inhomogenous Hawkes process

Here we fit a univariate inhomogenous marked Hawkes process where the conditional intensity function is given by

$$\lambda(t) = \mu(t) + \alpha \Sigma_{i:\tau_i<t}\text{exp}(-\beta * (t-\tau_i)) $$ The background $\mu(t)$ is time varying, rather than being constant.

The following example uses simulated data.

set.seed(1)
library(hawkesbow)
# Simulate a Hawkes process with mu = 1+sin(t), alpha=1, beta =2
times <- hawkesbow::hawkes(1000, fun=function(y) {1+0.5*sin(y)}, M=1.5, repr=0.5, family="exp", rate=2)$p

We will attempt to recover these parameter values, modelling the background as $ \mu(t) = A + Bsin(t)$. The background will be written as a function of $x$ and $y$, where $A = e^x$ and $B= logit(y) e^x$. This formulation ensures the background is never negative.

## The background function must take a single parameter and the time(s) at which it is evaluated
background <- function(params,times){
        A = exp(params[[1]])
        B = stats::plogis(params[[2]]) * A
        return(A + B*sin(times))
}

## The background_integral function must take a single parameter and the time at which it is evaluated
background_integral <- function(params,x){
        A = exp(params[[1]])
        B = stats::plogis(params[[2]]) * A
        return((A*x)-B*cos(x))
}
param = list(alpha = 0.5, beta = 1.5)
background_param = list(1,1)
fit <- fit_hawkes_cbf(times = times, parameters = param, background = background, background_integral = background_integral, background_parameters = background_param)

The estimated values of $A$ and $B$ respectively are

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

The estimated values of $\alpha$ and $\beta$ respectively are:

ab <- get_coefs(fit)[1:2,1]
ab

The sim_hawkes() function

args(sim_hawkes)

method = 1

sim <- sim_hawkes(mu = 2, alpha = 0.2, beta = 0.3, plot = TRUE)
head(sim)

method = 2

sim <- sim_hawkes(mu = 2, alpha = 0.2, beta = 0.3, plot = TRUE, method = 2)
head(sim)

cmjt/stelfi documentation built on Oct. 25, 2023, 2:53 p.m.