#' Self-exciting Hawkes process(es)
#'
#' Fit a Hawkes process using Template Model Builder (TMB). The function \code{fit_hawkes()} fits a
#' self-exciting Hawkes process with a constant background rate. Whereas, \code{fit_hawkes_cbf()} fits a Hawkes
#' processes with a user defined custom background function (non-homogeneous background rate). The function
#' \code{fit_mhawkes()} fits a multivariate Hawkes process that allows for between- and within-stream
#' self-excitement.
#'
#' @details A univariate Hawkes (Hawkes, AG. 1971) process is a self-exciting temporal point process
#' with conditional intensity function
#' \eqn{\lambda(t) = \mu + \alpha \Sigma_{i:\tau_i<t}e^{(-\beta * (t-\tau_i))}}. Here \eqn{\mu} is the constant
#' baseline rate, \eqn{\alpha} is the instantaneous increase in intensity after an event, and \eqn{\beta} is the
#' exponential decay in intensity. The term \eqn{\Sigma_{i:\tau_i<t} \cdots} describes the historic dependence
#' and the clustering density of the temporal point process, where the \eqn{\tau_i} are the events in
#' time occurring prior to time \eqn{t}. From this we can derive the following quantities 1) \eqn{\frac{\alpha}{\beta}}
#' is the branching ratio, it gives the average number of events triggered by an event, and
#' 2) \eqn{\frac{1}{\beta}} gives the rate of decay of the self-excitement.
#' Including mark information results in the conditional intensity
#' \eqn{\lambda(t; m(t)) = \mu + \alpha \Sigma_{i:\tau_i<t}m(\tau_i)e^{(-\beta * (t-\tau_i))}},
#' where \eqn{m(t)} is the temporal mark. This model can be fitted with \code{fit_hawkes()}.
#'
#' @references Hawkes, AG. (1971) Spectra of some self-exciting and mutually exciting point processes.
#' \emph{Biometrika}, \strong{58}: 83--90.
#'
#' @param times A vector of numeric observed time points.
#' @param parameters A named list of parameter starting values:
#' \itemize{
#' \item \code{mu}, base rate of the Hawkes process (\eqn{\mu}),
#' \item \code{alpha} (supplied if \code{model = 1}), intensity jump after an event occurrence (\eqn{\alpha}),
#' \item \code{a_par} (supplied if \code{model} = 2), logit scale for \code{alpha} (must be
#' set so that the intensity never goes negative and so that \code{alpha} <= \code{beta}), and
#' \item \code{beta}, exponential intensity decay (\eqn{\beta}).
#' }
#' @param model A numeric indicator specifying which model to fit:
#' \itemize{
#' \item \code{model = 1}, fits a Hawkes process with exponential decay (default);
#' \item \code{model = 2}, fits a Hawkes process with an \code{alpha} that can be negative.
#' }
#' @param marks Optional, a vector of numeric marks, defaults to 1 (i.e., no marks).
#' @param tmb_silent Logical, if \code{TRUE} (default) then
#' TMB inner optimisation tracing information will be printed.
#' @param optim_silent Logical, if \code{TRUE} (default) then for each iteration
#' \code{optim()} output will be printed.
#' @param ... Additional arguments to pass to \code{optim()}
#' @return A list containing components of the fitted model, see \code{TMB::MakeADFun}. Includes
#' \itemize{
#' \item \code{par}, a numeric vector of estimated parameter values;
#' \item \code{objective}, the objective function;
#' \item \code{gr}, the TMB calculated gradient function; and
#' \item \code{simulate}, (where available) a simulation function.
#' }
#' @examples
#' \donttest{
#' ### ********************** ###
#' ## A Hawkes model
#' ### ********************** ###
#' data(retweets_niwa, package = "stelfi")
#' times <- unique(sort(as.numeric(difftime(retweets_niwa, min(retweets_niwa), units = "mins"))))
#' params <- c(mu = 0.05, alpha = 0.05, beta = 0.1)
#' fit <- fit_hawkes(times = times, parameters = params)
#' get_coefs(fit)
#' ### ********************** ###
#' ## A Hawkes model with marks (ETAS-type)
#' ### ********************** ###
#' data("nz_earthquakes", package = "stelfi")
#' earthquakes <- nz_earthquakes[order(nz_earthquakes$origintime),]
#' earthquakes <- earthquakes[!duplicated(earthquakes$origintime), ]
#' times <- earthquakes$origintime
#' times <- as.numeric(difftime(times, min(times), units = "hours"))
#' marks <- earthquakes$magnitude
#' params <- c(mu = 0.05, alpha = 0.05, beta = 1)
#' fit <- fit_hawkes(times = times, parameters = params, marks = marks)
#' get_coefs(fit)
#' }
#' @seealso \code{\link{show_hawkes}}
#' @export
fit_hawkes <- function(times, parameters = list(), model = 1,
marks = c(rep(1, length(times))),
tmb_silent = TRUE,
optim_silent = TRUE, ...) {
## general error checks
for (i in 2:length(times)) {
if ((times[i] - times[i - 1]) < 1.e-10)
stop("times must be in ascending order with no simultaneous events")
}
if (length(marks) != length(times))
stop("marks must have same length as times")
if (min(marks) < 0)
stop("marks cannot be negative")
## parameters
alpha <- parameters[["alpha"]]
if (is.null(alpha)) {
alpha <- 0.5 * max(times) / length(times)
}
beta <- parameters[["beta"]]
if (is.null(beta)) {
beta <- 2 * alpha
}
mu <- parameters[["mu"]]
if (is.null(mu)) {
mu <- 0.5 / beta
}
if(model == 1) {
## error checks
if (alpha > (beta / mean(marks)))
stop("alpha must be smaller than or equal to beta divided by the mean of the marks")
if (alpha < 0)
stop("alpha must be non-negative")
## setup
obj <- TMB::MakeADFun(data = list(times = times, marks = marks, model_type = "hawkes"),
parameters = list(log_mu = log(mu),
logit_abratio = stats::qlogis(alpha / (beta / mean(marks))),
log_beta = log(beta)),
hessian = TRUE, DLL = "stelfi", silent = tmb_silent)
}else{
if (model == 2) {
## setup
obj <- TMB::MakeADFun(data = list(times = times, marks = marks, model_type = "neg_alpha_hawkes"),
parameters = list(log_mu = log(mu),
a_par = alpha,
log_beta = log(beta)),
hessian = TRUE, DLL = "stelfi", silent = tmb_silent)
}
}
trace <- if(optim_silent) 0 else 1
opt <- stats::optim(obj$par, obj$fn, obj$gr, control = list(trace = trace), ...)
obj$objective <- opt$value
return(obj)
}
#' @details An in-homogenous marked Hawkes process has conditional intensity function
#' \eqn{\lambda(t) = \mu(t) + \alpha \Sigma_{i:\tau_i<t}e^{(-\beta * (t-\tau_i))}}. Here, the
#' background rate, \eqn{\mu(t)}, varies in time. Such a model can be fitted
#' using \code{fit_hawkes_cbf()} where the parameters of the custom background function are estimated
#' before being passed to \code{TMB}.
#'
#' @param background A function taking one parameter and an
#' independent variable, returning a scalar.
#' @param background_integral The integral of \code{background}.
#' @param background_parameters The parameter(s)
#' for the background function \code{background}.
#' This could be a list of multiple values.
#' @param background_min A function taking one parameter and two points,
#' returns min of \code{background} between those points.
#' @examples
#' \donttest{
#' ### ********************** ###
#' ## A Hawkes process with a custom background function
#' ### ********************** ###
#' if(require("hawkesbow")) {
#' times <- hawkesbow::hawkes(1000, fun = function(y) {1 + 0.5*sin(y)},
#' M = 1.5, repr = 0.5, family = "exp", rate = 2)$p
#' ## 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)
#' get_coefs(fit)
#' }
#' }
#' @export
#' @rdname fit_hawkes
fit_hawkes_cbf <- function(times, parameters = list(),
model = 1,
marks = c(rep(1, length(times))),
background, background_integral, background_parameters,
background_min,
tmb_silent = TRUE, optim_silent = TRUE) {
## general error checks
for (i in 2:length(times)) {
if ((times[i] - times[i - 1]) < 1.e-10)
stop("times must be in ascending order with no simultaneous events")
}
if (length(marks) != length(times))
stop("marks must have same length as times")
if (min(marks) < 0)
stop("marks cannot be negative")
## beta parameter
beta <- parameters[["beta"]]
if (is.null(beta)) {
beta <- max(times) / length(times)
}
if (model == 1) {
## alpha parameter
alpha <- parameters[["alpha"]]
if (is.null(alpha)) {
alpha <- 0.5 * beta
}
## error checks
if (alpha > (beta / mean(marks)))
stop("alpha must be smaller than or equal to beta divided by the mean of the marks")
if (alpha < 0)
stop("alpha must be non-negative")
## Nested function to be passed into optim
optimize_background_one <- function(background_parameters, times, parameters,
marks, background, background_integral,
tmb_silent = TRUE, optim_silent = TRUE){
lambda <- background(background_parameters, times)
lambda_integral <- background_integral(background_parameters, tail(times, n = 1)) -
background_integral(background_parameters, 0)
alpha <- parameters[["alpha"]]
beta <- parameters[["beta"]]
obj <- TMB::MakeADFun(data = list(times = times, lambda = lambda,
lambda_integral = lambda_integral,
marks = marks, model_type = "custom_hawkes"),
parameters = list(logit_abratio = stats::qlogis(alpha / (beta / mean(marks))),
log_beta = log(beta)),
hessian = TRUE, DLL = "stelfi", silent = tmb_silent)
trace <- if(optim_silent) 0 else 1
opt <- stats::optim(obj$par, obj$fn, obj$gr, control = list(trace = 0))
return(opt$value)
}
opt <- stats::optim(par = background_parameters, fn = optimize_background_one,
gr = NULL, times, parameters, marks, background,
background_integral,
tmb_silent, optim_silent)
## Need to run again to extract alpha and beta
lambda <- background(opt$par, times)
lambda_integral <- background_integral(opt$par, tail(times, n = 1)) -
background_integral(opt$par, 0)
alpha <- parameters[["alpha"]]
beta <- parameters[["beta"]]
obj <- TMB::MakeADFun(data = list(times = times, lambda = lambda,
lambda_integral = lambda_integral,
marks = marks, model_type = "custom_hawkes"),
parameters = list(logit_abratio = stats::qlogis(alpha / beta),
log_beta = log(beta)),
hessian = TRUE, DLL = "stelfi", silent = tmb_silent)
}else{
if (model == 2){
## a_par parameter
a_par <- parameters[["a_par"]]
if (is.null(a_par)) {
a_par <- 0
}
## Nested function to be passed into optim
optimize_background_two <- function(background_parameters, times, parameters,
marks, background, background_integral,
background_min, tmb_silent = TRUE, optim_silent = TRUE) {
lambda <- background(background_parameters, times)
lambda_min <- numeric(length = length(times))
for (k in 1:(length(times) - 1)) {
lambda_min[k] <- background_min(background_parameters, times[k], times[k + 1])
}
lambda_min[length(times)] <- background(background_parameters, tail(times, n = 1))
if (min(lambda - lambda_min) < 0) stop("lambda_min incorrectly defined: lambda_min > lambda")
lambda_integral <- background_integral(background_parameters, tail(times, n = 1)) -
background_integral(background_parameters, 0)
a_par <- parameters[["a_par"]]
beta <- parameters[["beta"]]
obj <- TMB::MakeADFun(data = list(times = times, lambda = lambda,
lambda_min = lambda_min,
lambda_integral = lambda_integral,
marks = marks,
model_type = "neg_alpha_custom_hawkes"),
parameters = list(a_par = a_par,
log_beta = log(beta)),
hessian = TRUE, DLL = "stelfi", silent = tmb_silent)
trace <- if(optim_silent) 0 else 1
opt <- stats::optim(obj$par, obj$fn, obj$gr, control = list(trace = trace))
return(opt$value)
}
opt <- stats::optim(par = background_parameters, fn = optimize_background_two,
gr = NULL, times, parameters, marks, background,
background_integral,
background_min, tmb_silent, optim_silent)
## Need to run again to extract alpha and beta
lambda <- background(opt$par, times)
lambda_min <- numeric(length = length(times))
for (k in 1:(length(times) - 1)) {
lambda_min[k] <- background_min(opt$par, times[k], times[k + 1])
}
lambda_min[length(times)] <- background(opt$par, tail(times, n = 1))
lambda_integral <- background_integral(opt$par, tail(times, n = 1)) -
background_integral(opt$par, 0)
a_par <- parameters[["a_par"]]
beta <- parameters[["beta"]]
## setup
obj <- TMB::MakeADFun(data = list(times = times, lambda = lambda,
lambda_min = lambda_min,
lambda_integral = lambda_integral,
marks = marks,
model_type = "neg_alpha_custom_hawkes"),
parameters = list(a_par = a_par,
log_beta = log(beta)),
hessian = TRUE, DLL = "stelfi", silent = tmb_silent)
}
}
trace <- if(optim_silent) 0 else 1
opt2 <- stats::optim(obj$par, obj$fn, obj$gr, control = list(trace = trace))
obj$objective <- opt2$value
obj$background_parameters <- opt$par
return(obj)
}
#' @details A multivariate Hawkes process that allows for between- and within-stream self-excitement.
#' The conditional intensity for the \eqn{j^{th}} (\eqn{j = 1, ..., N}) stream is given by
#' \eqn{\lambda(t)^{j*} = \mu_j + \Sigma_{k = 1}^N\Sigma_{i:\tau_i<t} \alpha_{jk} e^{(-\beta_j * (t-\tau_i))}}, where
#' \eqn{j, k \in (1, ..., N)}. Here, \eqn{\alpha_{jk}} is the excitement caused by the
#' \eqn{k^{th}} stream on the \eqn{j^{th}}. Therefore, \eqn{\boldsymbol{\alpha}} is an \eqn{N x N}
#' matrix where the diagonals represent the within-stream excitement and the off-diagonals
#' represent the excitement between streams.
#'
#' @param stream A character vector specifying the stream ID of each observation in \code{times}
#' @param parameters A named list of parameter starting values:
#' \itemize{
#' \item \code{mu}, a vector of base rates for each stream of the multivariate Hawkes process,
#' \item \code{alpha}, a matrix of the between- and within-stream self-excitement: the diagonal
#' elements represent the within-stream excitement and the off-diagonals the excitement between streams,
#' \item \code{beta}, a vector of the exponential intensity decay for each stream of the multivariate
#' Hawkes process,
#' }
#' @examples
#' \donttest{
#' ### ********************** ###
#' ## A multivariate Hawkes model
#' ### ********************** ###
#' data(multi_hawkes, package = "stelfi")
#' fit <- 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)
#' }
#' @export
#' @rdname fit_hawkes
fit_mhawkes <- function(times, stream,
parameters = list(),
tmb_silent = TRUE,
optim_silent = TRUE, ...){
## general error checks
for (i in 2:length(times)) {
if ((times[i] - times[i - 1]) < 1.e-10)
stop("times must be in ascending order with no simultaneous events")
}
n_streams <- length(table(stream))
## parameters
alpha <- parameters[["alpha"]]
beta <- parameters[["beta"]]
mu <- parameters[["mu"]]
## error checks
if (n_streams != nrow(parameters$alpha) | n_streams != ncol(parameters$alpha))
stop("dimensions of alpha must match the number of streams")
if (sum(alpha > beta) > 0)
stop("alpha(s) must be smaller than or equal to beta")
if (sum(alpha < 0) > 0)
stop("alpha(s) must be non-negative")
events <- as.numeric(factor((stream))) - 1
events_per_stream <- as.numeric(table(stream)) - 1
obj <- TMB::MakeADFun(data = list(times = times, events = events,
N = n_streams,
EPS = events_per_stream,
model_type = "multi_hawkes"),
parameters = list(log_mu = log(mu),
logit_abratio = stats::qlogis(alpha/beta),
log_beta = log(beta)),
hessian = TRUE, DLL = "stelfi", silent = tmb_silent)
trace <- if(optim_silent) 0 else 1
opt <- stats::optim(obj$par, obj$fn, obj$gr, control = list(trace = trace), ...)
obj$objective <- opt$value
return(obj)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.