Reference

The public names, as listed in hawkes_engine.__all__. Anything else is private and may change without notice.

Estimation

hawkes_engine.fit(times, types, mu=None, alpha=None, beta=None, T=None, *, M=None, shared_decay=False, n_starts=1, init='moment', seed=None, se_method='hessian', ridge_report=True, zero_diagonal=False, n_bootstrap=200, auto_max_seconds=120.0, bootstrap_workers=None)

Fit the model by maximum likelihood; return the estimates and what is known about them.

The one constraint that matters: the reported standard errors are asymptotic and were measured to understate the sampling spread of the estimates. See hawkes_engine.trust_lookup() for the measured ratio in a regime, and the “Limits” page of the documentation.

Parameters:

• times (ArrayLike) – The event stream, as for hawkes_engine.loglik(): times sorted and non-negative, types integers in 0..M-1 (at most 256 types).

• types (ArrayLike) – The event stream, as for hawkes_engine.loglik(): times sorted and non-negative, types integers in 0..M-1 (at most 256 types).

• mu (ArrayLike | None) – Optional fixed values. A block supplied here is held constant and is not estimated; this is what makes a Poisson-limit test possible without a second entry point. None estimates that block.

• alpha (ArrayLike | None) – Optional fixed values. A block supplied here is held constant and is not estimated; this is what makes a Poisson-limit test possible without a second entry point. None estimates that block.

• beta (ArrayLike | None) – Optional fixed values. A block supplied here is held constant and is not estimated; this is what makes a Poisson-limit test possible without a second entry point. None estimates that block.

• T (float | None) – End of the observation window. Defaults to the last event time. Pass it explicitly when the window extends past the last event: the likelihood depends on T through the compensator.

• M (int | None) – Number of types. Inferred from mu or alpha when given, otherwise required.

• shared_decay (bool) – Estimate one decay per source rather than one per pair (ADR-002).

• n_starts (int) – Number of independent optimizer runs. The best log-likelihood wins.

• init (str) – "moment" (the only value in v0.1). See _moment_start().

• seed (int | None) – Seed for the multi-start generator and the bootstrap stream. None draws from the OS entropy pool, so a fit without a seed is not reproducible; with one it is bit-reproducible.

• se_method (str | None) – "auto" to let the calibration map choose (see hawkes_engine.se_policy and the se_source and se_note fields of the result); "hessian" for the observed information (the default, and the only method the calibration study measured); "bootstrap" for a parametric bootstrap that refits n_bootstrap times; or None to skip standard errors, in which case reliable is "unreliable" by construction. On Windows and macOS the bootstrap starts worker processes, so the calling script needs an if __name__ == "__main__": guard. Without one, Python 3.11 and later fail with an error that names it; Python 3.10 on Windows can hang instead (a CPython defect), so use the guard there. "auto" may choose the bootstrap, so a script that asks for "auto" needs the same guard.

• ridge_report (bool) – Accepted for compatibility. It has no effect: the reliability report always runs.

• zero_diagonal (bool) – Constrain alpha_ii = 0. Masked out of the free vector, not projected after the fact: a projected coordinate is still a coordinate.

• n_bootstrap (int) – Replications for se_method="bootstrap", and for the bootstrap that "auto" may choose.

• auto_max_seconds (float | None) – The ceiling on the projected cost of a bootstrap that "auto" would choose. When the projection exceeds it, "auto" keeps the Hessian interval, says so in se_note and warns. None removes the ceiling, which is a decision to let a fit call run for as long as the bootstrap takes. Only read when se_method="auto".

• bootstrap_workers (int | None) – Process count for the bootstrap pool. None uses every logical core, which is right for a single fit and wrong when many fits run at once: a study that parallelises over fits and leaves this at None starts one pool per fit and spends its time in process start-up. Set it to 1 there.

Returns:

The estimates, their standard errors where they exist, and a reliable verdict. reliable describes how well conditioned this fit is; it says nothing about bias or coverage, which are properties of the estimator and are reported by hawkes_engine.trust_lookup(). se_source and se_note say which method produced se and why; with se_method="auto" the auto field carries the whole decision.

Return type:

FitResult

Raises:

ValueError – If the event stream is unsorted, an event lies after T, a type index is not less than M, a supplied fixed block is outside its domain, se_method is not one of the three values above, or n_starts is below 1.

Notes

Convergence is judged on the gradient recomputed at the returned point, not on scipy’s status. A fit that stops on a search bound is reported as such via bound_hit rather than returned silently.

class hawkes_engine.FitResult(mu, alpha, beta, loglik, converged, status, nit, gradient_inf_norm, n_params, rho, rho_ci, se, cov, hessian, hessian_condition, reliable, notes=(), data=<factory>, bound_hit=False, bound_coordinates=(), theta=None, se_source='hessian', se_note='', auto=None)

A maximum-likelihood fit and everything known about its trustworthiness.

Parameters:

• mu (NDArray[float64])

• alpha (NDArray[float64])

• beta (NDArray[float64])

• loglik (float)

• converged (bool)

• status (int)

• nit (int)

• gradient_inf_norm (float)

• n_params (int)

• rho (float)

• rho_ci (tuple[float, float] | None)

• se (dict[str, NDArray[float64]] | None)

• cov (NDArray[float64] | None)

• hessian (NDArray[float64] | None)

• hessian_condition (float | None)

• reliable (str)

• notes (tuple[str, ...])

• data (dict[str, Any])

• bound_hit (bool)

• bound_coordinates (tuple[str, ...])

• theta (NDArray[float64] | None)

• se_source (str)

• se_note (str)

• auto (Any)

mu, alpha, beta

Natural-space estimates. alpha is (M, M) row = target, column = source (ADR-002). beta is (M, M) or (M,) when the decay is shared, in which case beta_ij = beta_j.

loglik

The log-likelihood at the estimate.

Type:

float

converged

The optimizer finished and the recomputed log-space gradient is within the convergence tolerance. status and nit are scipy’s.

Type:

bool

gradient_inf_norm

||d ell / d theta||_inf at the returned point, recomputed here rather than read from the optimizer’s last call.

Type:

float

rho

Spectral radius of the branching matrix Gamma = alpha / beta.

Type:

float

rho_ci

Delta-method 95% interval for rho, or None when the covariance is unavailable. It is a first-order interval and is not valid at a boundary; it is reported because rho is the headline number and an interval is more honest than a point even when it is approximate.

Type:

tuple[float, float] | None

se, cov

cov is the log-space covariance of the free parameters; se maps "mu", "alpha", "beta" to natural-space standard errors of the same shape as the estimates. Both are None when the information matrix is not invertible to the stated tolerance — never a number the library does not believe.

hessian

Observed information -d2 ell / d theta2 in natural space over the free parameters, in Layout order.

Type:

numpy._typing._array_like.NDArray[numpy.float64] | None

hessian_condition

cond(H) from the singular values, not the ratio of diagonal entries.

Type:

float | None

reliable

One of "reliable", "marginal", "unreliable", "unknown". This states whether the fit converged and is well conditioned. It does not state whether the estimator is unbiased or whether an interval covers at the nominal rate; those are properties of the estimator, not of this fit; the pre-registered study measured them (see hawkes_engine.trust_lookup()). The value comes from hawkes_engine.report(); "unknown" is a placeholder that fit replaces before returning.

Type:

str

notes

One line per reason a field is None, non-finite, or not to be trusted. Empty when there is nothing to say.

Type:

tuple[str, …]

data

The inputs and derived summaries: n_events, T, M, shared_decay, seed, zero_diagonal, n_starts, se_method, residual_gaps_summary.

Type:

dict[str, Any]

n_params

Number of free parameters.

Type:

int

status, nit

scipy’s termination status and iteration count for the winning start.

bound_hit, bound_coordinates

Whether a coordinate stopped on a search bound, and which. A fit that did is a boundary value, not a stationary point.

theta

The log-space optimum over the free parameters.

Type:

numpy._typing._array_like.NDArray[numpy.float64] | None

se_source

"hessian", "bootstrap" or "none": which method produced se. "none" means no interval was computed or the calibration map says the estimator is undefined in this regime. Read it before reading se.

Type:

str

se_note

One sentence saying why se_source has the value it has. Empty only for a result read from a 0.1.0 document, which did not record one.

Type:

str

auto

The full AutoDecision record when se_method="auto", else None. It names the map cell, the distance, the q it read, the threshold and the projected cost, so the choice can be audited from the result alone.

Type:

Any

classmethod from_dict(payload)

Rebuild from to_dict(). Rejects a schema it does not know.

Parameters:

payload (dict[str, Any])

Return type:

FitResult

to_dict()

A JSON-serialisable mapping. Stable: see SCHEMA_VERSION.

Return type:

dict[str, Any]

Likelihood

hawkes_engine.loglik(times, types, mu, alpha, beta, T, *, grad=False, shared_decay=False)

Evaluate the log-likelihood of a multivariate Hawkes process.

Implements docs/spec/mathematics.md section 3:

\[\ell(\theta) = \sum_i \left[ \sum_k \log \lambda_i(t^i_k) - \int_0^T \lambda_i(t)\,dt \right]\]

Parameters:

• times (ArrayLike) – Event times, shape (N,), non-decreasing, seconds from the start of the session. Convert epoch timestamps to a session-relative origin before calling (ADR-003).

• types (ArrayLike) – Type index of each event, shape (N,), in 0..M-1.

• mu (ArrayLike) – Baseline intensities, shape (M,), non-negative.

• alpha (ArrayLike) – Excitation block, shape (M, M), non-negative. Row is the target type, column the source (ADR-002).

• beta (ArrayLike) – Decay rates, shape (M, M) for a distinct decay per pair, or (M,) when shared_decay is true, in which case beta_ij = beta_j. All entries strictly positive.

• T (float) – End of the observation window. Must be at least the last event time.

• grad (bool) – Also compute the gradient with respect to the natural-space parameters.

• shared_decay (bool) – Use the source-only decay parameterization, O(N*M) instead of O(N*M*M).

Returns:

value and, when requested, dmu, dalpha, dbeta.

Return type:

LoglikResult

Raises:

ValueError – If the events are not sorted or lie outside [0, T], a parameter is outside its domain, or a shape disagrees with M.

Notes

The gradient is in natural space. The optimizer’s log-space gradient is natural * parameter entrywise; that conversion is applied once, in the objective wrapper, not here (ADR-002).

Two events at the same timestamp do not excite each other: the intensity sums over strictly earlier events.

class hawkes_engine.LoglikResult

Return value of loglik().

value is always present. The gradient blocks are present only when grad=True: dmu has shape (M,), dalpha and the pair dbeta have shape (M, M), and the shared-decay dbeta has shape (M,).

When the intensity is not strictly positive at some event, value is -inf and every gradient entry is nan. Nothing is clipped: a fit that reaches that region has found something wrong, and it should be able to see it.

Simulation

hawkes_engine.simulate(mu, alpha, beta, T, *, method='thinning', seed, burn_in=0.0, t0=0.0)

Simulate events on [t0, T); return (times, types).

Parameters:

• mu (ArrayLike) – Baseline intensities, shape (M,), non-negative.

• alpha (ArrayLike) – Excitation block, shape (M, M). Row is the target type, column the source (ADR-002).

• beta (ArrayLike) – Decay rates, shape (M, M), strictly positive.

• T (float) – End of the simulated window. Must be finite and greater than t0.

• method (Literal['thinning', 'cluster']) – "thinning" for Ogata’s modified thinning, "cluster" for the branching representation. Both sample the same law; neither is a fallback for the other.

• seed (int) – Seed of the stream, in [0, 2**64). The output is a function of the seed alone, so the same arguments reproduce the same path on every platform.

• burn_in (float) – Time to run before t0 and discard. Zero starts from an empty history, which leaves the intensity below its stationary level early in the window.

• t0 (float) – Start of the returned window. Events before it are simulated and discarded when burn_in > 0.

Returns:

times (float64, sorted, strictly inside the window) and types (uint8), the same length. Both are writable and independent of the arguments.

Return type:

tuple of numpy.ndarray

Raises:

ValueError – When the spectral radius of alpha / beta is 1 or more, the window is malformed, or the parameters are outside their domain.

hawkes_engine.simulate_many(mu, alpha, beta, T, *, n_paths, method='thinning', seed, burn_in=0.0, t0=0.0)

Simulate n_paths independent paths.

Parameters:

• mu (ArrayLike) – As for simulate(). The parameter arrays are converted once and reused, so the per-path cost is the simulation.

• alpha (ArrayLike) – As for simulate(). The parameter arrays are converted once and reused, so the per-path cost is the simulation.

• beta (ArrayLike) – As for simulate(). The parameter arrays are converted once and reused, so the per-path cost is the simulation.

• T (float) – As for simulate(). The parameter arrays are converted once and reused, so the per-path cost is the simulation.

• method (Literal['thinning', 'cluster']) – As for simulate(). The parameter arrays are converted once and reused, so the per-path cost is the simulation.

• seed (int) – As for simulate(). The parameter arrays are converted once and reused, so the per-path cost is the simulation.

• burn_in (float) – As for simulate(). The parameter arrays are converted once and reused, so the per-path cost is the simulation.

• t0 (float) – As for simulate(). The parameter arrays are converted once and reused, so the per-path cost is the simulation.

• n_paths (int) – Number of paths. A negative value raises ValueError; zero returns an empty list.

Returns:

One (times, types) pair per path.

Return type:

list of tuple of numpy.ndarray

Notes

Path k uses the seed seed_stream(seed, k), so the result is a function of (seed, k) alone and does not depend on how many threads the caller uses. The paths are produced serially here; a caller who splits the index range across threads gets the same per-path seeds.

Diagnostics

hawkes_engine.time_rescale(times, types, mu, alpha, beta)

Rescaled gaps for each type, in event order.

Parameters:

• times (ArrayLike) – Event times, shape (N,), non-decreasing.

• types (ArrayLike) – Type index of each event, shape (N,), in 0..M-1.

• mu (ArrayLike) – Parameters to rescale under, with the shapes hawkes_engine.loglik() documents. M is the length of mu.

• alpha (ArrayLike) – Parameters to rescale under, with the shapes hawkes_engine.loglik() documents. M is the length of mu.

• beta (ArrayLike) – Parameters to rescale under, with the shapes hawkes_engine.loglik() documents. M is the length of mu.

Returns:

One float64 array per type. Under the true parameters each is independent Exp(1) (mathematics.md section 8), which is what ks_test() and ljung_box() test.

Return type:

list of numpy.ndarray

Notes

The first gap of each type runs from time 0. The stretch after a type’s last event is not a gap and is dropped. A type with no events yields an empty array, not an error.

Raises ValueError if the events are unsorted or a type index is not less than M.

hawkes_engine.ks_test(gaps)

One-sample KS test of each type’s gaps against Exp(1).

Parameters:

gaps (list[NDArray[float64]]) – Rescaled gaps, one array per type, as returned by time_rescale().

Returns:

statistic and pvalue, one entry per element of gaps.

Return type:

TestResult

Notes

Per type, so an empty type gets nan rather than 1.0: there is nothing to test, and a pass would read as evidence. Use it together with ljung_box(); the marginal law can be right while the dependence is wrong.

The p-values assume the parameters were fixed before seeing the data. When the gaps come from parameters fitted to the same events, the nominal p-values are not calibrated, and a pass is absence of evidence against the model, not evidence for it.

hawkes_engine.qq_data(gaps)

(theoretical, sample) quantile pairs against Exp(1), per type.

Parameters:

gaps (list[NDArray[float64]]) – Rescaled gaps, one array per type, as returned by time_rescale().

Returns:

Per type, the theoretical quantiles and the sorted sample. Plot them with matplotlib.pyplot.plot; this module does not import pyplot, so a caller who never draws a plot never loads it.

Return type:

list of tuple of numpy.ndarray

Notes

Theoretical quantiles use the plotting positions (k - 1/2) / n, which keep every point inside (0, 1) so the logarithm is finite.

hawkes_engine.report(model, n_events=None)

Assess conditioning from this fit alone.

Parameters:

• model (FitResult) – A fitted FitResult. Only fields the fit already computed are read; nothing is refitted.

• n_events (int | None) – Event count to compare against the parameter count. None reads model.data["n_events"], which is what the fit recorded.

Returns:

A verdict of "reliable", "marginal" or "unreliable", the evidence, and one note per reason.

Return type:

ReliabilityReport

Notes

Reads the Hessian condition number, whether a coordinate stopped on a search bound, the branching ratio, the alpha-beta correlation implied by the covariance, and the event count per parameter. It reports whether the fit converged and is well-conditioned, not whether it is unbiased.

Trust

hawkes_engine.trust_lookup(N, rho_hat, beta, T, resolution=0.0, *, M=2, rows=None)

The measured cell nearest to a fitted model’s regime.

The one constraint that matters: the answer describes the nearest cell of one fixed parameter set per cell, not the query’s own parameters, and it is given however far that cell is. Check distance.

Parameters:

• N (float) – Number of events the model was fitted to.

• rho_hat (float) – The fitted branching ratio.

• beta (float) – The smallest fitted decay rate, the decay scale of the slowest kernel.

• T (float) – The observation window.

• resolution (float) – Timestamp resolution in the same units as T; 0 for exact times. The measured range of beta * resolution is 0 to 0.13.

• M (int) – Number of event types. Only 2 (and 4, the satellite) were measured.

• rows (list[dict[str, str]] | None) – Map rows to search; None loads the shipped map. For tests.

Returns:

status is "measured", "unknown" or "estimator undefined in this regime". Coverage and q are given only when measured, each as (estimate, Monte Carlo SE) per block, with the distance to the cell they were measured in.

Return type:

MapAnswer

class hawkes_engine.MapAnswer(status, cell, distance, coverage, q, convergence_rate, no_se_rate, study_id, note)

What the map can say about a regime, and how far away it measured it.

Parameters:

• status (str)

• cell (str | None)

• distance (float | None)

• coverage (dict[str, tuple[float, float]] | None)

• q (dict[str, tuple[float, float]] | None)

• convergence_rate (float | None)

• no_se_rate (float | None)

• study_id (str)

• note (str)

status

"measured", "unknown" or "estimator undefined in this regime".

Type:

str

cell

Key of the nearest measured grid cell, or None when status is "unknown".

Type:

str | None

distance

Distance to that cell in grid-step units over (log10 N, rho, r, log10 beta_min*T); 0 is a grid point, 1 is about one grid step. The lookup does not refuse a large distance, so read it.

Type:

float | None

coverage

Per block ("mu", "alpha", "beta", "all"), the measured coverage of nominal-95 percent intervals as (estimate, Monte Carlo SE). A missing standard error counts as a miss. None unless measured.

Type:

dict[str, tuple[float, float]] | None

q

Per block, (mean reported SE / empirical SD, bootstrap SE). Below 1 means the reported standard errors were smaller than the spread of the estimates. None unless measured.

Type:

dict[str, tuple[float, float]] | None

convergence_rate, no_se_rate

Share of replications in the cell that converged, and that returned no standard error. None when status is "unknown".

study_id, note

The study that produced the map, and a sentence on how to read the answer.

Stability

hawkes_engine.spectral_radius(mu, alpha, beta)

rho(Gamma) for Gamma = alpha / beta: the branching ratio.

Parameters:

• mu (ArrayLike) – Baseline intensities, shape (M,). Unused by the computation but part of the parameter set, so it is required and validated here rather than silently defaulted.

• alpha (ArrayLike) – Excitation and decay blocks, shape (M, M), with beta strictly positive.

• beta (ArrayLike) – Excitation and decay blocks, shape (M, M), with beta strictly positive.

Returns:

The Perron root of the entrywise non-negative matrix Gamma, certified to a bracket width below 1e-13.

Return type:

float

Notes

A value at or above 1 means the process is supercritical: no stationary distribution exists, and simulation refuses the parameters. The function itself reports the value for any parameters, because knowing how far past criticality a set is is often the question.

hawkes_engine.stationary_mean(mu, alpha, beta)

(I - Gamma)^{-1} mu: the mean intensity of the stationary process.

Parameters:

• mu (ArrayLike) – As for spectral_radius().

• alpha (ArrayLike) – As for spectral_radius().

• beta (ArrayLike) – As for spectral_radius().

Returns:

Shape (M,), float64. A new array, writable, with no reference to the arguments.

Return type:

numpy.ndarray

Raises:

ValueError – When rho(Gamma) >= 1: the inverse does not exist as a positive quantity, and the answer would be a meaningless sign pattern rather than a mean.