Multivariate Hawkes Processes

• Multivariate Hawkes processes are point process models whose intensity functions depend on past events using self- and cross-excitation kernels.
• The process structure features a baseline vector and a kernel matrix that encode directed causal interactions and stability through a branching representation.
• Modern extensions incorporate marked, spatial, inhibition, and neural formulations to capture complex event dependencies and enhance scalability.

Multivariate Hawkes processes are multivariate point processes whose predictable conditional intensities depend on past events across multiple coordinates through self- and cross-excitation kernels. In the classical linear formulation, a $d$-dimensional process $N(t)=(N_1(t),\dots,N_d(t))$ has intensities

$$\lambda_i(t)=\nu_i+\sum_{j=1}^d\int_{(0,t)} h_{ij}(t-s)\,dN_j(s),$$

so the matrix of kernels encodes directed interaction, stability, and, in the mutually exciting case, a Poisson cluster interpretation (Boly et al., 2023). Contemporary work extends this template to marked, discrete-time, spatio-temporal, nonlinear, quadratic, and neural formulations, including models with common drivers, simultaneous multi-coordinate events, and explicit inhibition (Lotz, 2024, Bielecki et al., 2020, Chukwuemeka et al., 27 Feb 2026).

1. Formal structure and canonical representations

A multivariate Hawkes process is usually specified by a baseline vector and a matrix of causal kernels. In the standard unmarked case, the integrated kernel matrix

$B_{ij}=\int_0^\infty \phi_{ij}(u)\,du$

is the branching or adjacency matrix. Its entries quantify the expected number of type-$i$ offspring generated by a type-$j$ event, and, under the point-process Granger-causality interpretation, the directed edge $j\to i$ is present if and only if the corresponding interaction coefficient is nonzero (Lotz, 2024). This support-based interpretation is mirrored in the nonparametric graphical formulation, where $j$ does not Granger-cause $i$ exactly when the link function $h_{ij}(\cdot)$ is identically zero (Eichler et al., 2016).

Marked formulations enlarge this structure by allowing event-specific attributes to modulate excitation. A marked multivariate Hawkes process may be written as

$N(t)=(N_1(t),\dots,N_d(t))$0

or, in separable form, $N(t)=(N_1(t),\dots,N_d(t))$1, so marks scale the excitation while the temporal kernel controls decay (Lotz, 2024). A complementary perspective is to represent a univariate marked Hawkes process as a multivariate unmarked process by partitioning the mark space into bins $N(t)=(N_1(t),\dots,N_d(t))$2 and treating each bin as a component. The induced intensity takes the form

$N(t)=(N_1(t),\dots,N_d(t))$3

which yields an $N(t)=(N_1(t),\dots,N_d(t))$4 approximation for a broad class of marked Hawkes processes and preserves stationarity when the target process is stationary (Davis et al., 2024).

Spatial structure can be added by replacing event times $N(t)=(N_1(t),\dots,N_d(t))$5 with $N(t)=(N_1(t),\dots,N_d(t))$6 and considering

$N(t)=(N_1(t),\dots,N_d(t))$7

This formulation admits separable and nonseparable kernels, and it underlies both classical spatio-temporal Hawkes models and recent neural generalizations (Chukwuemeka et al., 27 Feb 2026).

2. Stability, branching structure, and asymptotic regimes

For linear mutually exciting multivariate Hawkes processes with nonnegative integrable kernels, the standard stationarity condition is spectral: if

$N(t)=(N_1(t),\dots,N_d(t))$8

then the stationary mean intensity exists and is

$N(t)=(N_1(t),\dots,N_d(t))$9

Under the same condition, the process admits a Poisson cluster decomposition in which immigrants arrive as independent Poisson processes and each cluster is a subcritical multitype Galton–Watson branching process with offspring mean matrix $$\lambda_i(t)=\nu_i+\sum_{j=1}^d\int_{(0,t)} h_{ij}(t-s)\,dN_j(s),$$0 (Boly et al., 2023).

This classical subcritical regime supports sharp dependence and limit theory. If the kernels have finite $$\lambda_i(t)=\nu_i+\sum_{j=1}^d\int_{(0,t)} h_{ij}(t-s)\,dN_j(s),$$1-moments, the stationary multivariate linear Hawkes process is strongly mixing with polynomial decay,

$$\lambda_i(t)=\nu_i+\sum_{j=1}^d\int_{(0,t)} h_{ij}(t-s)\,dN_j(s),$$2

and a functional central limit theorem holds for linear functionals and for the centered counting process. The long-run covariance matrix of the vector limit is

$$\lambda_i(t)=\nu_i+\sum_{j=1}^d\int_{(0,t)} h_{ij}(t-s)\,dN_j(s),$$3

which is also the Bartlett spectrum at frequency zero (Boly et al., 2023).

Near criticality, the relevant asymptotics are different. For asymptotically critical sequences of multivariate Hawkes processes, suitable time-space scaling yields weak convergence to a multidimensional stochastic Volterra equation with measure kernel,

$$\lambda_i(t)=\nu_i+\sum_{j=1}^d\int_{(0,t)} h_{ij}(t-s)\,dN_j(s),$$4

where the potential measure $$\lambda_i(t)=\nu_i+\sum_{j=1}^d\int_{(0,t)} h_{ij}(t-s)\,dN_j(s),$$5 is associated with a matrix-valued extended Bernstein function determined by the scaled kernel asymptotics (Xu, 2024). In the continuous-limit case, the martingale part becomes a time-changed Brownian motion,

$$\lambda_i(t)=\nu_i+\sum_{j=1}^d\int_{(0,t)} h_{ij}(t-s)\,dN_j(s),$$6

placing asymptotically critical Hawkes systems within the class of affine Volterra processes (Xu, 2024).

Large-network limits on graphs lead to a further notion of stability. On inhomogeneous random graphs, the macroscopic nonlinear limit is the neural-field-type equation

$$\lambda_i(t)=\nu_i+\sum_{j=1}^d\int_{(0,t)} h_{ij}(t-s)\,dN_j(s),$$7

while in the linear case the long-time regime is governed by the spectral radius $$\lambda_i(t)=\nu_i+\sum_{j=1}^d\int_{(0,t)} h_{ij}(t-s)\,dN_j(s),$$8 of the graphon operator $$\lambda_i(t)=\nu_i+\sum_{j=1}^d\int_{(0,t)} h_{ij}(t-s)\,dN_j(s),$$9. The subcritical threshold is $B_{ij}=\int_0^\infty \phi_{ij}(u)\,du$0, and the supercritical regime $B_{ij}=\int_0^\infty \phi_{ij}(u)\,du$1 yields growth of the macroscopic intensity norm (Agathe-Nerine, 2021).

3. Causality, graphs, marks, and common drivers

The causal interpretation of multivariate Hawkes processes is unusually direct. In the classical linear framework, the support of the kernel matrix coincides with the Granger-causality graph: $B_{ij}=\int_0^\infty \phi_{ij}(u)\,du$2 if and only if $B_{ij}=\int_0^\infty \phi_{ij}(u)\,du$3 (Eichler et al., 2016). In the parametric marked setting, the same role is played by the adjacency coefficients $B_{ij}=\int_0^\infty \phi_{ij}(u)\,du$4, where

$B_{ij}=\int_0^\infty \phi_{ij}(u)\,du$5

Testing whether a specified group of $B_{ij}=\int_0^\infty \phi_{ij}(u)\,du$6 vanishes is therefore a formal test of causal influence (Lotz, 2024).

This support-based viewpoint motivates both estimation and hypothesis testing. When many independent short trajectories are available, the aggregation device

$B_{ij}=\int_0^\infty \phi_{ij}(u)\,du$7

turns the sample into a Hawkes process with large baseline, enabling constrained maximum likelihood under positivity and zero restrictions on the tested edges. Because the null lies on the boundary $B_{ij}=\int_0^\infty \phi_{ij}(u)\,du$8, the likelihood-ratio statistic converges to a chi-bar-squared mixture rather than a standard $B_{ij}=\int_0^\infty \phi_{ij}(u)\,du$9 law (Lotz, 2024).

Marks can either be treated explicitly or absorbed into a multivariate representation. Partition-based multivariate representations show that continuous or discrete mark dependence can be approximated by cross-excitation across mark bins, with identifiable parameters under mild conditions and stationarity preserved when the target marked process is stationary (Davis et al., 2024). This suggests that “multivariate” need not only mean multiple event types in the raw data; it may also be an inferential representation of heterogeneous marks.

A related decomposition separates global and local sources of dependence. In relational event networks with covariates, one may write

$i$0

where $i$1 and $i$2 are common drivers shared across actors, while $i$3 and $i$4 are local parameters. This specification is designed to distinguish synchronous activity induced by a common driver from genuine network excitation (Kreiss et al., 4 Apr 2025).

4. Estimation, scalable learning, and exact sampling

Likelihood-based inference is canonical but computationally demanding. For linear Hawkes processes with general kernels, both the log-likelihood and the least-squares error, as well as their gradients, have quadratic complexity in the number of observed events under naive evaluation (Cartea et al., 2021). A key scalable alternative is the least-squares decomposition based on the primitives

$i$5

which allows unbiased adaptive stratified Monte Carlo estimates of the gradient. The resulting ASLSD algorithm reduces the per-iteration complexity to

$i$6

for kernel families with closed-form $i$7 and $i$8 (Cartea et al., 2021).

Dimension scaling motivates structural reductions. The Low-Rank Hawkes Process parameterizes

$i$9

so the $j$0 kernels are replaced by $j$1 latent interactions and a $j$2 nonnegative projection matrix. With an exponential basis, the resulting learning complexity is

$j$3

rather than $j$4 (Lemonnier et al., 2016). A different route exploits sparsity of individual cascades rather than low rank of the interaction matrix. In the “Lazy MHP” construction, the cross-excitation coefficients factor as $j$5 for $j$6, which allows exact likelihoods and gradients whose cost depends on the number of active entities in each sequence rather than on the ambient network size (Nickel et al., 2020).

Bayesian scalable inference has also been developed. For exponential-kernel MHPs, stochastic gradient EM, stochastic gradient variational inference, and stochastic gradient Langevin Monte Carlo can all be implemented with contiguous time-window minibatches. A boundary-corrected compensator approximation preserves conjugacy-like structure in SGEM and SGVI while reducing approximation errors associated with right-boundary effects in short windows (Jiang et al., 2023).

Simulation and sampling are equally varied. Ogata-style thinning can be adapted to space-time and marked settings, but exact stationary sampling is possible in the mutually exciting multivariate case through a perfect sampling construction based on immigrant–offspring clusters, exponential tilting of cluster birth times, and accept–reject correction. The algorithm generates i.i.d. stationary sample paths without transient bias, and its expected complexity has an explicit convex form in the tilting parameters (Chen et al., 2020).

5. Beyond mutual excitation: inhibition, simultaneity, quadratic feedback, and neural state dynamics

The classical linear model is not the full theory. A generalized multivariate Hawkes process may allow inhibitory kernels and nonlinear link functions through

$j$7

where $j$8 may take negative values and $j$9 is only required to be Lipschitz. Existence and uniqueness are obtained under the contraction condition $j\to i$0, with

$j\to i$1

thereby extending high-dimensional Hawkes theory beyond purely excitatory linear kernels (Chen et al., 2017).

Another extension targets explicit simultaneous multi-coordinate events. Generalized multivariate Hawkes processes define marks in

$j\to i$2

so a single event time can carry active marks in several coordinates at once. This construction yields common event times with positive probability and separates idiosyncratic excitation from explicitly joint excitation channels (Bielecki et al., 2020).

Quadratic Hawkes processes alter the feedback mechanism itself. In the multivariate quadratic extension,

$j\to i$3

so intensities depend not only on past event counts but also on quadratic interactions in past returns. With the ZHawkes decomposition, the quadratic kernel splits into a time-diagonal Hawkes term and a rank-one trend contribution, allowing endogenous co-jumps and heavy-tailed volatility distributions (Aubrun et al., 2022).

Neural parameterizations replace explicit kernels with latent state dynamics. The multivariate spatio-temporal neural Hawkes process specifies

$j\to i$4

with a continuous-time LSTM state whose memory cell decays according to

$j\to i$5

This parameterization learns excitation and inhibition while preserving nonnegativity through the softplus link, and it avoids predefined triggering kernels (Chukwuemeka et al., 27 Feb 2026).

6. Discrete-time formulations, applications, and recurrent limitations

Not all event data are naturally continuous-time. In discrete Hawkes models, counts in each bin are conditionally Poisson, and the multivariate intensity becomes a lagged count recursion. For the marked multivariate discrete formulation used in incident monitoring,

$j\to i$6

which is more appropriate when timestamps are rounded and simultaneous events are common (Brisley et al., 2023).

Applications in the supplied literature span finance, epidemiology, neural spike trains, social information diffusion, earthquakes, order-book dynamics, terrorism, and incident monitoring. Formal testing of Hawkes causal graphs has been illustrated on retail online auctions and German intraday power prices (Lotz, 2024). A marked discrete multivariate Hawkes process with alarm marks was used for violent incidents in a forensic psychiatric hospital, where the marked multivariate model achieved the best predictive log-likelihood among the compared alternatives (Brisley et al., 2023). In multivariate spatio-temporal neural Hawkes modeling, an application to terrorism data from Pakistan showed distinct temporal and spatial interaction patterns across four groups and highlighted learned inhibitory “cool-down” phases after bursts (Chukwuemeka et al., 27 Feb 2026). A general non-Markovian multivariate Hawkes framework with time-dependent baseline and arbitrary matrix-valued excitation kernel has also been proposed for spatio-temporal epidemics, together with multi-temporal Laplace transforms, moment formulas, and two-time covariance decompositions (Vieille et al., 21 Jul 2025).

Several recurrent limitations appear across these formulations. One is misspecification: collapsing intrinsically spatio-temporal data into purely temporal streams can severely distort recovered intensities (Chukwuemeka et al., 27 Feb 2026). Another is identifiability or parameter proliferation: multivariate representations of marked processes introduce $j\to i$7 parameters for exponential-kernel representations, and larger $j\to i$8 increases approximation flexibility at the cost of estimation difficulty (Davis et al., 2024). High-dimensional causal testing involves chi-bar-squared weights that are generally intractable beyond simple cases (Lotz, 2024). Classical mutually exciting models also do not cover inhibition or simultaneous multi-coordinate events without explicit generalization (Chen et al., 2017, Bielecki et al., 2020).

A persistent methodological theme is that likelihood alone is not sufficient as an evaluation criterion. Parametric, neural, and spatio-temporal studies alike emphasize direct analysis of fitted intensities, covariance structure, or recovered interaction shape, because high likelihood can coexist with implausible dynamics under misspecification or excessive flexibility (Chukwuemeka et al., 27 Feb 2026, Boly et al., 2023). This suggests that multivariate Hawkes processes are best understood not as a single model class but as a family of interacting point-process constructions whose common core is history-dependent multivariate intensity, and whose modern variants differ mainly in how they encode interaction, nonlinearity, marks, space, and computational tractability.