--- title: "Get started with emhawkes" author: "Kyungsub Lee" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteEncoding{UTF-8} %\VignetteIndexEntry{Get started with emhawkes} %\VignetteEngine{knitr::rmarkdown} editor_options: markdown: wrap: sentence --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) knitr::opts_chunk$set(eval = requireNamespace("maxLik", quietly = TRUE)) ``` This vignette introduces model specification, simulation, estimation and a first residual check. It is included in the R package and also published on the website. For more models, see [Extended examples](https://ksublee.github.io/emhawkes/articles/extended_example.html). ## Basic Hawkes model ### Univariate Hawkes process ```{r} library(emhawkes) ``` This section shows how to specify, simulate, and fit a univariate Hawkes model. First, create an `hspec` object to define the model. The S4 class `hspec` includes slots such as `mu`, `alpha`, `beta`, `dimens`, `rmark`, and `impact`. In a univariate model, `mu`, `alpha`, and `beta` can be supplied as numeric values, which are converted to matrices. The following example defines an unmarked univariate Hawkes model. ```{r} mu1 <- 0.3; alpha1 <- 1.2; beta1 <- 1.5 hspec1 <- new("hspec", mu = mu1, alpha = alpha1, beta = beta1) show(hspec1) ``` The `hsim` function simulates a path from an `hspec` object. The `size` argument specifies the number of rows, including the initial state. The optional arguments `lambda_component0` and `N0` specify the initial intensity components and event counts, respectively. The intensity of the basic univariate Hawkes model is $$ \lambda(t) = \mu + \int_{-\infty}^t \alpha e^{-\beta (t-s)} d N(s) = \mu + \lambda_c(0) e^{-\beta t} + \int_0^t \alpha e^{-\beta (t-s)} d N(s) $$ where `lambda_component0` denotes $$ \lambda_c(0) = \int_{-\infty}^0 \alpha e^{\beta s} d N(s). $$ The default `lambda_component0 = "auto"` uses unconditional component means for supported stable constant numeric linear models, and zero components otherwise. The selected rule is reported in a message; zero initialization recommends supplying components explicitly. Use `"mean"` to require supported finite means, `"zero"` to start from baseline intensities, or a scalar/matrix to specify the components directly. Mean initialization is deterministic and does not sample a stationary initial state. For a sufficiently long simulation of a stable model, the effect of the initial intensity diminishes over time. The default initial value of the counting process, `N0`, is zero. ```{r, warning=FALSE} set.seed(1107) res1 <- hsim(hspec1, size = 1000) summary(res1) ``` The `hsim` function returns an S3 object of class `hreal` with the following components. - `hspec` is the model specification. - `inter_arrival` contains the times between consecutive events. - `arrival` is the cumulative sum of `inter_arrival`. - `type` identifies the counting process associated with each event: type $i$ corresponds to $N_i$. - `mark` contains the mark associated with each event. - `N` contains the cumulative event counts for each type. - `Nc` contains the cumulative marks for each type. - `lambda` contains the intensity immediately before each event. - `rambda` contains the intensity immediately after each event. - `lambda_component` and `rambda_component` contain the intensity components $\lambda_{ij}$ immediately before and after each event, respectively. With the default initial values, `inter_arrival`, `type`, `mark`, `N`, and `Nc` begin with zero. The `summary()` function displays the first 20 rows of the realized path by default. The `print()` function also displays the model specification and intensity components. In this univariate model, `lambda == mu + lambda_component`. ```{r} # The first and third columns are equal. cbind(res1$lambda[1:5], res1$lambda_component[1:5], mu1 + res1$lambda_component[1:5]) ``` For all rows except the first, `rambda` equals `lambda + alpha` in this model. ```{r} # The second and third columns are equal except in the initial row. cbind(res1$lambda[1:5], res1$rambda[1:5], res1$lambda[1:5] + alpha1) ``` The following calculation verifies the exponential decay between events. ```{r} # By definition, the following two are equal: res1$lambda[2:6] mu1 + (res1$rambda[1:5] - mu1) * exp(-beta1 * res1$inter_arrival[2:6]) ``` Use `logLik` to evaluate the log-likelihood from the model specification and observed inter-arrival times. ```{r, warning=FALSE} logLik(hspec1, inter_arrival = res1$inter_arrival) ``` Use `hfit` for maximum likelihood estimation. The `hspec0` object supplies the initial parameter values for optimization. For this univariate model, `inter_arrival` is the only required data argument. If the initial intensity component is known, pass it as `lambda_component0`; otherwise, the function determines it internally. The default optimization method is BFGS. ```{r, warning=FALSE} # Initial parameter values for numerical optimization. mu0 <- 0.5; alpha0 <- 1.0; beta0 <- 1.8 hspec0 <- new("hspec", mu = mu0, alpha = alpha0, beta = beta0) # Supply the initial values through hspec0. mle <- hfit(hspec0, inter_arrival = res1$inter_arrival) summary(mle) ``` Use `fitted_hspec()` to reconstruct the estimated model and `fit_initial_state()` to reuse its initial-state settings. The fit remains a `maxLik` object. `coef()`, `summary()`, and `vcov()` provide the usual estimation results. ```{r fit-reconstruction, message=FALSE} estimated <- fitted_hspec(mle) initial <- fit_initial_state(mle) inferred <- do.call(logLik, c(list(object = estimated, inter_arrival = res1$inter_arrival, infer = TRUE), initial)) c(fitted = as.numeric(logLik(mle)), restored = inferred$loglikelihood) head(inferred$lambda) head(inferred$res_process1) ``` ### A first residual check For an ordinary Hawkes model, integrated-intensity residuals should follow an exponential distribution with rate 1 when the model is adequate. Compare their quantiles with the reference distribution. ```{r residual-check, fig.height=4, fig.width=4} p <- ppoints(100) q <- quantile(inferred$res_process1, probs = p) plot(qexp(p), q, xlab = "Theoretical quantiles", ylab = "Sample quantiles") abline(0, 1, col = "blue", lty = 2) ``` By default, `hfit()` uses the final Hessian and standard errors from `maxLik`. For saved fits, new simulations, fuller residual diagnostics and optional `final_hessian = "stable"` checks, see [Fitting and inference](https://ksublee.github.io/emhawkes/articles/fitting-and-inference.html). A stable Hessian alone does not establish reliable confidence intervals. ### Bivariate Hawkes model The intensity process of a basic bivariate Hawkes model is defined by $$ \lambda_1(t) = \mu_1 + \int_{-\infty}^t \alpha_{11} e^{-\beta_{11}(t-s)} d N_1(s) + \int_{-\infty}^t \alpha_{12} e^{-\beta_{12}(t-s)} d N_2(s), $$ $$ \lambda_2(t) = \mu_2 + \int_{-\infty}^t \alpha_{21} e^{-\beta_{21}(t-s)} d N_1(s) + \int_{-\infty}^t \alpha_{22} e^{-\beta_{22}(t-s)} d N_2(s). $$ For this bivariate model, `mu` is a 2-by-1 matrix, and `alpha` and `beta` are 2-by-2 matrices. $$ \mu = \begin{bmatrix} \mu_1 \\ \mu_2 \end{bmatrix}, \quad \alpha = \begin{bmatrix} \alpha_{11} & \alpha_{12} \\ \alpha_{21} & \alpha_{22} \end{bmatrix}, \quad \beta = \begin{bmatrix} \beta_{11} & \beta_{12} \\ \beta_{21} & \beta_{22} \end{bmatrix} $$ `rmark` generates marks during simulation and can be omitted for unmarked models. The 2-by-2 matrix `lambda_component0` specifies the initial values of `lambda11`, `lambda12`, `lambda21`, and `lambda22`. The intensity processes are represented by $$ \lambda_1(t) = \mu_1 + \lambda_{11}(t) + \lambda_{12}(t), $$ $$ \lambda_2(t) = \mu_2 + \lambda_{21}(t) + \lambda_{22}(t). $$ The terms $\lambda_{ij}$ are called intensity components, and `lambda_component0` specifies their initial values $\lambda_{ij}(0)$. If `lambda_component0` is omitted, the function determines these values internally. ```{r} mu2 <- matrix(c(0.2), nrow = 2) alpha2 <- matrix(c(0.5, 0.9, 0.9, 0.5), nrow = 2, byrow = TRUE) beta2 <- matrix(c(2.25, 2.25, 2.25, 2.25), nrow = 2, byrow = TRUE) hspec2 <- new("hspec", mu=mu2, alpha=alpha2, beta=beta2) print(hspec2) ``` Use `hsim` to simulate the model. ```{r} set.seed(1107) res2 <- hsim(hspec2, size=1000) summary(res2) ``` In multivariate models, `type` identifies the counting process associated with each event. ```{r} # A bivariate model has two event types: 1 and 2. res2$type[1:10] ``` In multivariate models, the column names of `N` are `N1`, `N2`, `N3`, and so on. ```{r} res2$N[1:3, ] ``` Similarly, the column names of `lambda` are `lambda1`, `lambda2`, `lambda3`, and so on. ```{r} res2$lambda[1:3, ] ``` In this bivariate model, the columns of `lambda_component` are `lambda11`, `lambda12`, `lambda21`, and `lambda22`. ```{r} res2$lambda_component[1:3, ] ``` By definition, the following two expressions are equivalent: ```{r} mu2[1] + rowSums(res2$lambda_component[1:5, c("lambda11", "lambda12")]) res2$lambda[1:5, "lambda1"] ``` Extract the observed `inter_arrival` and `type` vectors from the simulation result. Both are required to fit a bivariate model. ```{r} inter_arrival2 <- res2$inter_arrival type2 <- res2$type ``` Use `logLik` to evaluate the log-likelihood. ```{r} logLik(hspec2, inter_arrival = inter_arrival2, type = type2) ``` Use `hfit` to estimate the parameters from `inter_arrival` and `type`. The values of `mu`, `alpha`, and `beta` in `hspec0` provide starting points for optimization. Here, we use the simulation parameters by setting `hspec0 <- hspec2`. For observed data, choose initial values because the true parameters are unknown. By default, `reduced = TRUE` groups equal starting values within each parameter matrix. Here, the two entries of `mu` share one parameter, the diagonal and off-diagonal entries of `alpha` give two parameters, and all entries of `beta` share one parameter. The fit therefore estimates four parameters. `hfit` prints a short explanation once before optimization. ```{r, warning=FALSE} hspec0 <- hspec2 mle <- hfit(hspec0, inter_arrival = inter_arrival2, type = type2) summary(mle) ``` ```{r} coef(mle) ``` ```{r} miscTools::stdEr(mle) ``` The restored matrices retain the shared parameters from this fit. The diagonal entries of `alpha` remain equal, as do its off-diagonal entries. ```{r shared-matrices} estimated2 <- fitted_hspec(mle) estimated2@alpha estimated2@beta ``` Set `reduced = FALSE` to estimate all ten matrix entries separately, even when their starting values are equal. Use `suppressMessages(hfit(...))` to hide the startup explanation. For function components, parameters with the same name are shared under either setting. A zero in a numeric matrix is an estimated starting value; use a component function to fix selected entries at zero. ### Observed data and initial state The first row represents the initial state and is not an event. For event times measured from the observation start, prepare the input as follows: ```{r, eval=FALSE} inter_arrival <- c(0, diff(c(0, event_times))) type <- c(0, event_types) mark <- c(0, event_marks) ``` Omit `type` for a univariate model and `mark` for unit marks. The vectors must have the same length; event types range from 1 to the model dimension. Event order is preserved, including events at the same time. The likelihood currently ends at the last supplied event. `N0` and `Nc0` specify initial event counts and mark totals. Pass `Nc0` by name in `hfit`, and use the same initial state when comparing simulation and likelihood results. `lambda_component0` can be a scalar or a matrix matching `alpha` and `beta`. If supplied, `N` and `Nc` must have one row per input row and one column per dimension, and must agree with the events and initial counts. Their rows contain the cumulative values after each event. ## Choosing the optimizer `hfit()` uses the optional maxLik package by default. To use R's optimization functions directly, select `optimizer = "stats"`; likelihood engine selection and model specification are unchanged. This route also works without maxLik. The earlier fitting examples require maxLik and are not evaluated when it is unavailable. ```{r stats-optimizer, eval=TRUE} library(emhawkes) poisson <- new("hspec", mu = 0.5, alpha = 0, beta = 1) fit_stats <- hfit(poisson, c(0, 0.2, 0.3, 0.5), optimizer = "stats", lambda_component0 = 0, fixed = c(FALSE, TRUE, TRUE), control = list(maxit = 500, reltol = 1e-12)) summary(fit_stats) fit_diagnostics(fit_stats)$computation$optimizer ``` The stats route supports BFGS, Nelder-Mead, CG and L-BFGS-B. Linear inequality constraints use `constrOptim()` and require an interior starting point; equality constraints are not supported. Use `control$maxit` rather than maxLik's `iterlim`. It defaults to stable final curvature; boundary estimates have no Wald covariance. Use `fitted_hspec()` and `fit_initial_state()` with either provider. Changing providers does not guarantee identical estimates or cure convergence and identification problems. ## Further reading - [Fitting and inference](https://ksublee.github.io/emhawkes/articles/fitting-and-inference.html): parameter settings, initial states, fitted models, residual diagnostics and standard errors. - [Extended examples](https://ksublee.github.io/emhawkes/articles/extended_example.html): marked, state-dependent, multi-kernel and flexible residual models. - [Computation](https://ksublee.github.io/emhawkes/articles/computation.html): likelihood engines, analytic gradients, Jacobians and performance. - [Reference](https://ksublee.github.io/emhawkes/reference/index.html): function arguments and return values.