Loading required package: stelfi
mark | jump |
---|---|
1 | 0.4 |
2 | 0.8 |
3 | 1.2 |
4 | 1.6 |
The conditional intensity function for the marked Hawkes model implemented in stelfi
is given by
Below data are simulated using the emhawkes
package (Lee 2023) where the marks (a vector that scale the jump sizes, starting at 0) are integer values stelfi
). The parameter values of the conditional intensity are
Loading required package: stelfi
mark | jump |
---|---|
1 | 0.4 |
2 | 0.8 |
3 | 1.2 |
4 | 1.6 |
To fit the model in stelfi
the fit_hawkes()
function is used and the additional optional argument marks
supplied.
sv <- c(mu = 1.3, alpha = 0.4, beta = 1.5)
fit <- fit_hawkes(times = res$arrival, parameters = sv, marks = res$mark)
get_coefs(fit)
Estimate Std. Error
mu 1.1741551 0.25195698
alpha 0.1008830 0.06935281
beta 0.8849185 0.50179668
Note the estimated coefficient stelfi
equates to emhawkes
.
## benchmark emhawkes
emhawkes::hfit(h1, inter_arrival = res$inter_arrival, mark = res$mark) |>
summary()
--------------------------------------------
Maximum Likelihood estimation
BFGS maximization, 43 iterations
Return code 0: successful convergence
Log-Likelihood: -51.99802
3 free parameters
Estimates:
Estimate Std. error t value Pr(> t)
mu1 1.1106 0.2888 3.845 0.00012 ***
alpha1 0.2787 0.1903 1.464 0.14309
beta1 0.9277 0.5693 1.629 0.10321
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
--------------------------------------------
The fitted model and diagnostic plots are plotted using show_hawkes()
and show_hawkes_GOF()
.
TODO
Below the negative log-likelihood of a univariate marked Hawkes process used by stelfi
is written using R
syntax and then RTMB
(Kristensen 2024) is used for parameter estimation. See here for and overview of RTMB
. This section is for demonstration only, feel free to modify the function as desired. Note that this can be used to fit an unmarked model by setting the vector of marks to be
library(RTMB)
univariate_marked_hawkes <- function(params){
getAll(data, params)
mu <- exp(log_mu)
beta <- exp(log_beta)
alpha <- exp(logit_abratio) / (1 + exp(logit_abratio)) * (beta/mean(marks))
n <- length(times)
last <- times[n]
nll <- 0
A <- advector(numeric(n))
for(i in 2:n){
A[i] <- sum(exp(-beta * (times[i] - times[i - 1])) * (marks[i - 1] + A[i - 1]))
}
term_3vec <- log(mu + alpha * A)
nll <- (mu * last) + ((alpha/beta) * (sum(marks) - marks[n] - A[n])) - sum(term_3vec)
ADREPORT(mu)
ADREPORT(alpha)
ADREPORT(beta)
return(nll)
}
data <- list(times = res$arrival, marks = res$mark)
params <- list(log_mu = log(1.3), logit_abratio = 0.6, log_beta = log(1.5))
obj <- MakeADFun(univariate_marked_hawkes, params, silent = TRUE)
opt <- nlminb(obj$par, obj$fn, obj$gr)
summary(sdreport(obj), "report")
Estimate Std. Error
mu 1.1742764 0.25196080
alpha 0.1008268 0.06934413
beta 0.8846088 0.50173890