Concepts

The notation is that of the specification, which is normative and holds the derivations. This page states what the code computes and what it costs.

The model

There are \(M\) event types. The intensity of type \(i\) at time \(t\) is

\[\lambda_i(t) = \mu_i + \sum_{j=1}^{M} \sum_{t^j_m < t} \alpha_{ij}\, e^{-\beta_{ij}(t - t^j_m)}.\]

\(\mu_i\) is a baseline rate, \(\alpha_{ij}\) the jump a type-\(j\) event adds to the intensity of type \(i\), and \(\beta_{ij}\) the rate at which that jump decays. Only exponential kernels and non-negative \(\alpha\) are supported. The sum is over strictly earlier events, so two events at the same timestamp do not excite each other.

Index convention. Row is the target, column is the source: alpha[i, j] acts from type j onto type i. This is the transpose of some references and is fixed by ADR-002. A transposed alpha runs without error and describes a different process.

Branching ratio and stationarity

The branching matrix is \(\Gamma_{ij} = \alpha_{ij} / \beta_{ij}\), the expected number of type-\(i\) children of one type-\(j\) event. The process is stationary when the spectral radius \(\rho(\Gamma) < 1\), and its mean intensity is then \((I - \Gamma)^{-1}\mu\). spectral_radius and stationary_mean compute both from parameters alone. simulate refuses parameters with \(\rho \ge 1\). fit does not constrain \(\rho\) to be below 1 and reports whatever it finds.

The log-likelihood and the one-pass recursion

For events observed on \([0, T]\),

\[\ell(\theta) = \sum_{i=1}^{M}\left[\sum_k \log \lambda_i(t^i_k) - \int_0^T \lambda_i(t)\,dt\right].\]

Evaluated directly this is quadratic in the number of events \(N\). Exponential kernels allow a recursion: keep the \(M^2\) states \(S_{ij}(t) = \sum_{t^j_m < t} e^{-\beta_{ij}(t - t^j_m)}\), decay all of them by \(e^{-\beta_{ij}\Delta}\) between consecutive events, and add 1 to \(S_{id}\) for every \(i\) when an event of type \(d\) arrives. The likelihood and its gradient come out of a single pass over the sorted events.

Cost. One evaluation is \(O(N M^2)\) time and \(O(M^2)\) extra memory for pair decay, and \(O(NM)\) when the decay depends only on the source (shared_decay=True). The measured cost per event, fit times and memory are on Performance. The reference implementation the compiled code is checked against is the naive quadratic one.

Two simulators

simulate has two algorithms, method="thinning" (Ogata’s modified thinning) and method="cluster" (the branching representation: Poisson immigrants, each event producing Poisson(\(\Gamma_{ij}\)) children after Exp(\(\beta_{ij}\)) delays). They are written separately and share only the random-number generator and the argument checks. That is the reason there are two: each is tested against the analytic stationary mean and against the other, so a defect in one is not hidden by the same defect in the other. Both sample the same law; neither is a fallback.

Both start from an empty history unless burn_in is given, which leaves the intensity below its stationary level early in the window. The calibration study discarded ten relaxation times, \(10 / (\beta_{\min}(1 - \rho))\).

Time rescaling

If the parameters are right, the integrated intensity between consecutive events of each type is an independent Exp(1) variable (Papangelou). time_rescale returns those gaps per type. The first gap of a type runs from time 0, and the stretch after its last event is dropped. ks_test checks the marginal law and hawkes_engine.diagnostics.ljung_box checks independence; the first can pass while the second fails.

Time origin

Event times are floats. Epoch seconds near \(1.7 \times 10^9\) leave about \(2 \times 10^{-7}\) s of resolution, so subtract the session start before calling (ADR-003).