Hierarchical Hawkes Process Models

• Hierarchical Hawkes process models extend traditional self-exciting point processes by imposing multiple layers on baselines, kernels, and latent variables to capture complex event dynamics.
• They integrate Bayesian priors, network factorization, and nonstationary segmentation techniques to improve parameter estimation and infer latent structures in multivariate systems.
• These models find applications in seismology, finance, social networks, and epidemics, while also confronting challenges such as computational complexity and identifiability issues.

Searching arXiv for recent and foundational papers on hierarchical and related Hawkes-process formulations. Searching arXiv for "hierarchical Hawkes process multivariate Bayesian segmentation network". Searching arXiv for relevant Hawkes-process surveys and hierarchy-related models. arXiv search query: Hawkes processes survey hierarchy multivariate Bayesian segmentation network. A hierarchical Hawkes process denotes a Hawkes-process construction in which the conditional intensity retains the self-exciting Hawkes form while multi-level structure is imposed on baselines, excitation kernels, latent variables, marks, or interaction networks; in a distinct but related usage, it denotes a hierarchical multi-resolution segmentation of a single nonstationary Hawkes process rather than a unified generative hierarchy (Laub et al., 2024, Zhou et al., 2019). All such models inherit the defining Hawkes property that each arrival increases the rate of future arrivals for some time, and they therefore rest on the same point-process, compensator, branching, and likelihood machinery as standard univariate and multivariate Hawkes processes (Laub et al., 2015).

1. Formal Hawkes foundations

A point process on $[0,\infty)$ may be specified by event times $\{T_1 \le T_2 \le \dots\}$ or by the counting process $N(t)$, the number of points in $(0,t]$. Its key local characterization is the conditional intensity

$$\lambda^*(t) = \lim_{h\downarrow 0} \frac{\mathbb{E}[N(t+h)-N(t)\mid H(t)]}{h},$$

with compensator

$\Lambda(t)=\int_0^t \lambda^*(s)\,ds.$

When it exists, $\lambda^*(t)$ determines all finite-dimensional distributions of the point process (Laub et al., 2015).

The linear univariate Hawkes process has intensity

$$\lambda^*(t)=\mu+\int_0^t \phi(t-u)\,dN(u) =\mu+\sum_{t_i<t}\phi(t-t_i),$$

where $\mu>0$ is the baseline intensity and $\phi(\cdot)\ge 0$ is the excitation kernel. Standard kernels include the exponential form

$\{T_1 \le T_2 \le \dots\}$0

and the power-law form

$\{T_1 \le T_2 \le \dots\}$1

The exponential kernel yields Markov dynamics for the pair $\{T_1 \le T_2 \le \dots\}$2, while the power-law kernel is associated with Omori-type aftershock decay (Laub et al., 2015).

The multivariate Hawkes process generalizes self-excitation to mutual excitation. For $\{T_1 \le T_2 \le \dots\}$3 counting processes $\{T_1 \le T_2 \le \dots\}$4,

$\{T_1 \le T_2 \le \dots\}$5

where $\{T_1 \le T_2 \le \dots\}$6 governs self-excitation and $\{T_1 \le T_2 \le \dots\}$7 with $\{T_1 \le T_2 \le \dots\}$8 governs cross-excitation. This matrix of kernels is the basic substrate for hierarchical structure, because it can be pooled, factorized, regularized, or linked through latent variables (Laub et al., 2015).

The immigration–birth representation interprets immigrants as a homogeneous Poisson process of rate $\{T_1 \le T_2 \le \dots\}$9, with each event independently generating offspring according to an inhomogeneous Poisson process with intensity $N(t)$0 for $N(t)$1. The scalar branching ratio is

$N(t)$2

and for the exponential kernel $N(t)$3. The univariate process is stable when $N(t)$4, with long-run mean intensity

$N(t)$5

In the multivariate case, defining $N(t)$6 gives a branching matrix $N(t)$7, and the stability condition becomes $N(t)$8, where $N(t)$9 is the spectral radius (Laub et al., 2015).

2. Meanings of hierarchy in Hawkes models

The phrase “hierarchical Hawkes process” is not used explicitly in the major surveys, but the relevant constructions fall into four recurrent classes: hierarchical Bayesian Hawkes, process-on-process hierarchies, network or factor hierarchies, and marked or spatio-temporal hierarchies (Laub et al., 2024). In all of them, the lower level consists of point-process event times generated by Hawkes intensities, while an upper level supplies shared parameters, latent regimes, latent factors, marks, or exogenous drivers.

Hierarchical Bayesian Hawkes models place higher-level prior distributions on baselines, kernel amplitudes, decay rates, or network structure across multiple related processes. Process-on-process hierarchies instead let one stochastic process modulate another, as in renewal immigration, shot-noise Cox baselines, or stochastic excitation amplitudes. Network or factor hierarchies impose structured dependence on the excitation matrix itself, through sparsity, low rank, clusters, or latent graphs. Marked and spatio-temporal hierarchies treat magnitude, type, or location as higher-level variables that organize offspring production and propagation (Laub et al., 2024, Lima, 2020).

A distinct usage arises in nonstationary Hawkes modeling. There, “hierarchical” refers to a nested sequence of temporal segmentations for a single process. The process is assumed piecewise stationary, and one constructs a coarse-to-fine tree of stationary regimes by ranking adjacent-sector discrepancies in estimated cumulants. This yields a hierarchical temporal decomposition, but not an explicit hierarchical generative model with linked priors across levels (Zhou et al., 2019).

3. Hierarchical Bayesian and factorized multivariate constructions

A generic hierarchical Bayesian Hawkes specification begins with unit-level intensities. For unit $(0,t]$0 and component $(0,t]$1,

$(0,t]$2

and with exponential kernels,

$(0,t]$3

A standard hierarchical prior pattern is

$(0,t]$4

$(0,t]$5

with hyperpriors on $(0,t]$6 and $(0,t]$7. This induces partial pooling and sharing of statistical strength across individuals and components (Laub et al., 2024).

Group-structured multivariate Hawkes models refine this by mapping indices to groups. One example is

$(0,t]$8

where $(0,t]$9 maps component $$\lambda^*(t) = \lim_{h\downarrow 0} \frac{\mathbb{E}[N(t+h)-N(t)\mid H(t)]}{h},$$0 to a group and $$\lambda^*(t) = \lim_{h\downarrow 0} \frac{\mathbb{E}[N(t+h)-N(t)\mid H(t)]}{h},$$1 is a truncated normal on $$\lambda^*(t) = \lim_{h\downarrow 0} \frac{\mathbb{E}[N(t+h)-N(t)\mid H(t)]}{h},$$2. Baseline intensities may similarly share common priors across related entities. Because stability requires $$\lambda^*(t) = \lim_{h\downarrow 0} \frac{\mathbb{E}[N(t+h)-N(t)\mid H(t)]}{h},$$3, the stable region becomes part of the prior support or parameterization, for instance through transformations or rejection of unstable proposals (Laub et al., 2015).

A complementary route separates structural connectivity from temporal response. One factorization writes

$$\lambda^*(t) = \lim_{h\downarrow 0} \frac{\mathbb{E}[N(t+h)-N(t)\mid H(t)]}{h},$$4

where $$\lambda^*(t) = \lim_{h\downarrow 0} \frac{\mathbb{E}[N(t+h)-N(t)\mid H(t)]}{h},$$5 is the impact matrix and $$\lambda^*(t) = \lim_{h\downarrow 0} \frac{\mathbb{E}[N(t+h)-N(t)\mid H(t)]}{h},$$6 contains kernel shapes. Low-rank models make the hierarchy explicit. In the Scalable Low-Rank Hawkes Process,

$$\lambda^*(t) = \lim_{h\downarrow 0} \frac{\mathbb{E}[N(t+h)-N(t)\mid H(t)]}{h},$$7

$$\lambda^*(t) = \lim_{h\downarrow 0} \frac{\mathbb{E}[N(t+h)-N(t)\mid H(t)]}{h},$$8

so that

$$\lambda^*(t) = \lim_{h\downarrow 0} \frac{\mathbb{E}[N(t+h)-N(t)\mid H(t)]}{h},$$9

Here the latent layer is an $\Lambda(t)=\int_0^t \lambda^*(s)\,ds.$0-dimensional Hawkes network and the observed $\Lambda(t)=\int_0^t \lambda^*(s)\,ds.$1 processes inherit intensities through the loading matrix $\Lambda(t)=\int_0^t \lambda^*(s)\,ds.$2 (Lima, 2020).

Mixture-based hierarchy appears in the Dirichlet–Hawkes process, where Hawkes components are drawn from a Dirichlet process prior and event sequences are modeled as superpositions of those components. This is a direct population-level distribution over Hawkes parameters, with cluster-level Hawkes dynamics nested beneath it (Lima, 2020).

4. Marked, spatio-temporal, and process-on-process hierarchies

Marked Hawkes processes extend event times by marks $\Lambda(t)=\int_0^t \lambda^*(s)\,ds.$3. The conditional intensity on $\Lambda(t)=\int_0^t \lambda^*(s)\,ds.$4 can be written

$\Lambda(t)=\int_0^t \lambda^*(s)\,ds.$5

with ground intensity

$\Lambda(t)=\int_0^t \lambda^*(s)\,ds.$6

This form is intrinsically hierarchical: the mark distribution defines an upper layer, while the ground Hawkes process determines event times. Marks may encode magnitude, type, size, location, cluster identity, or community. In this setting, the mark-dependent weight $\Lambda(t)=\int_0^t \lambda^*(s)\,ds.$7 modulates the offspring intensity induced by each event (Laub et al., 2024).

The ETAS model is a canonical marked hierarchical Hawkes model. Magnitudes follow the Gutenberg–Richter law

$\Lambda(t)=\int_0^t \lambda^*(s)\,ds.$8

while expected offspring productivity obeys Utsu’s law $\Lambda(t)=\int_0^t \lambda^*(s)\,ds.$9. The corresponding marked intensity is

$\lambda^*(t)$0

Large-magnitude events therefore generate more offspring, yielding a direct hierarchy from magnitude to offspring rate (Laub et al., 2024).

Spatio-temporal Hawkes models treat location as a mark or structured coordinate. In space–time ETAS,

$\lambda^*(t)$1

where $\lambda^*(t)$2 is a spatial kernel. This introduces a natural spatial hierarchy such as country, region, subregion, and local cluster, with magnitude and space jointly shaping local offspring intensity (Laub et al., 2024).

Process-on-process hierarchy is also explicit in several modern generalizations. Renewal Hawkes processes replace Poisson immigration by a renewal process. Dynamic contagion processes superpose external shot-noise excitation and internal Hawkes self-excitation:

$\lambda^*(t)$3

Stochastic excitation-kernel models let the excitation amplitudes $\lambda^*(t)$4 evolve as an SDE. In each case, one stochastic layer drives or modulates another (Laub et al., 2024).

5. Likelihood, inference, and computational constraints

For a point process observed on $\lambda^*(t)$5 with event times $\lambda^*(t)$6, the likelihood is

$\lambda^*(t)$7

and the log-likelihood is

$\lambda^*(t)$8

For a multivariate Hawkes process with intensities $\lambda^*(t)$9 and event times $$\lambda^*(t)=\mu+\int_0^t \phi(t-u)\,dN(u) =\mu+\sum_{t_i<t}\phi(t-t_i),$$0,

$$\lambda^*(t)=\mu+\int_0^t \phi(t-u)\,dN(u) =\mu+\sum_{t_i<t}\phi(t-t_i),$$1

With exponential kernels, the relevant inner sums can be updated recursively, so likelihood and derivatives can be computed in $$\lambda^*(t)=\mu+\int_0^t \phi(t-u)\,dN(u) =\mu+\sum_{t_i<t}\phi(t-t_i),$$2 time per evaluation (Laub et al., 2015).

The branching representation supplies latent immigrant/offspring structure and is central to inference. For a univariate Hawkes process, the probability that event $$\lambda^*(t)=\mu+\int_0^t \phi(t-u)\,dN(u) =\mu+\sum_{t_i<t}\phi(t-t_i),$$3 is an immigrant is

$$\lambda^*(t)=\mu+\int_0^t \phi(t-u)\,dN(u) =\mu+\sum_{t_i<t}\phi(t-t_i),$$4

while the probability that it is a child of a previous event $$\lambda^*(t)=\mu+\int_0^t \phi(t-u)\,dN(u) =\mu+\sum_{t_i<t}\phi(t-t_i),$$5 is

$$\lambda^*(t)=\mu+\int_0^t \phi(t-u)\,dN(u) =\mu+\sum_{t_i<t}\phi(t-t_i),$$6

EM algorithms use these parentage probabilities in the E-step and weighted parameter updates in the M-step. In hierarchical Bayesian settings, the same latent structure appears inside Gibbs sampling, Metropolis–Hastings, hybrid MCMC, RJMCMC for nonparametric kernels, and variational inference (Lima, 2020).

Hierarchical models sharpen classical Hawkes constraints rather than replacing them. Stability conditions such as $$\lambda^*(t)=\mu+\int_0^t \phi(t-u)\,dN(u) =\mu+\sum_{t_i<t}\phi(t-t_i),$$7 or $$\lambda^*(t)=\mu+\int_0^t \phi(t-u)\,dN(u) =\mu+\sum_{t_i<t}\phi(t-t_i),$$8 must be enforced by prior design or parameterization. High-dimensional multivariate models also create substantial statistical pressure: one financial example cited in the survey has $$\lambda^*(t)=\mu+\int_0^t \phi(t-u)\,dN(u) =\mu+\sum_{t_i<t}\phi(t-t_i),$$9 parameters for $\mu>0$0 assets, which raises overfitting concerns. Model selection criteria such as AIC and BIC, generalized method of moments, low-rank structure, shrinkage, and group regularization are therefore natural complements to hierarchical modeling (Laub et al., 2015).

Goodness-of-fit remains compensator-based. If $\mu>0$1 is the fitted intensity, transformed event times

$\mu>0$2

should form a unit-rate Poisson process under a correct model. In hierarchical settings this naturally extends to posterior predictive checks across units or groups (Laub et al., 2015).

6. Multi-resolution temporal hierarchy in nonstationary Hawkes models

The stationary Hawkes process assumes a constant baseline and a kernel depending only on lag. In the nonstationary setting studied by the multi-resolution segmentation literature, the baseline $\mu>0$3 is piecewise constant in time and the kernel becomes time-varying, $\mu>0$4, but is assumed locally stationary inside segments. Formally, there exists a partition of $\mu>0$5 into segments such that within each segment the process is a stationary Hawkes process with its own $\mu>0$6 and $\mu>0$7 (Zhou et al., 2019).

The segmentation machinery is built from first- and second-order cumulants. For a stationary Hawkes process,

$\mu>0$8

and the normalized conditional covariance is

$\mu>0$9

The kernel and $\phi(\cdot)\ge 0$0 satisfy the Wiener–Hopf equation

$\phi(\cdot)\ge 0$1

Within each finest-grid sector $\phi(\cdot)\ge 0$2, $\phi(\cdot)\ge 0$3 is estimated by a histogram

$\phi(\cdot)\ge 0$4

and adjacent sectors are compared through the normalized mean squared error

$\phi(\cdot)\ge 0$5

The algorithm sorts these discrepancy scores and, for a desired resolution $\phi(\cdot)\ge 0$6, selects the top $\phi(\cdot)\ge 0$7 boundaries as cut points (Zhou et al., 2019).

As $\phi(\cdot)\ge 0$8 increases from $\phi(\cdot)\ge 0$9 to $\{T_1 \le T_2 \le \dots\}$00, the segmentations are nested: $\{T_1 \le T_2 \le \dots\}$01 gives no cuts, while $\{T_1 \le T_2 \le \dots\}$02 cuts at every sector boundary. The resulting object is a coarse-to-fine temporal hierarchy. Technically, this is a global ranking of local discrepancy scores rather than recursive binary splitting, but coarser segments are unions of finer segments, so the output is naturally interpreted as a hierarchy of temporal regimes (Zhou et al., 2019).

A Gaussian-process refinement, GP-MRS, smooths the empirical $\{T_1 \le T_2 \le \dots\}$03 using

$\{T_1 \le T_2 \le \dots\}$04

with squared exponential kernel

$\{T_1 \le T_2 \le \dots\}$05

and posterior mean

$\{T_1 \le T_2 \le \dots\}$06

Plain cumulant estimation has complexity $\{T_1 \le T_2 \le \dots\}$07 and the boundary scoring adds $\{T_1 \le T_2 \le \dots\}$08, absorbed by $\{T_1 \le T_2 \le \dots\}$09. GP-MRS adds $\{T_1 \le T_2 \le \dots\}$10 per sector, yielding total complexity $\{T_1 \le T_2 \le \dots\}$11 (Zhou et al., 2019).

On the New York City vehicle-collision data, GP-MRS with $\{T_1 \le T_2 \le \dots\}$12, $\{T_1 \le T_2 \le \dots\}$13, $\{T_1 \le T_2 \le \dots\}$14, $\{T_1 \le T_2 \le \dots\}$15, and $\{T_1 \le T_2 \le \dots\}$16 revealed explicit nested regimes. For weekdays, the low-resolution structure at $\{T_1 \le T_2 \le \dots\}$17 separated busy and non-busy periods, while at $\{T_1 \le T_2 \le \dots\}$18 the cut positions were $\{T_1 \le T_2 \le \dots\}$19, $\{T_1 \le T_2 \le \dots\}$20, $\{T_1 \le T_2 \le \dots\}$21, and $\{T_1 \le T_2 \le \dots\}$22, with baseline intensities

$\{T_1 \le T_2 \le \dots\}$23

The nonstationary nonparametric Hawkes model achieved the lowest test negative log-likelihood $\{T_1 \le T_2 \le \dots\}$24 on held-out data, outperforming both stationary parametric and stationary nonparametric alternatives (Zhou et al., 2019).

7. Rank-structured network Hawkes models, applications, and limitations

A concrete hierarchical Hawkes model appears in network interaction data. For a directed social network $\{T_1 \le T_2 \le \dots\}$25, one may define a point process for each ordered pair $\{T_1 \le T_2 \le \dots\}$26, where events are wins of $\{T_1 \le T_2 \le \dots\}$27 over $\{T_1 \le T_2 \le \dots\}$28. The latent hierarchy is represented by continuous ranks

$\{T_1 \le T_2 \le \dots\}$29

with larger $\{T_1 \le T_2 \le \dots\}$30 indicating higher dominance. In the cohort Hawkes process, the dyad intensity is

$\{T_1 \le T_2 \le \dots\}$31

where

$\{T_1 \le T_2 \le \dots\}$32

This produces directional hierarchy at the kernel level: if $\{T_1 \le T_2 \le \dots\}$33, then $\{T_1 \le T_2 \le \dots\}$34 (Ward et al., 2020).

The most expressive model in that sequence is the cohort Markov-Modulated Hawkes process (C-MMHP). Each dyad has a latent two-state continuous-time Markov chain $\{T_1 \le T_2 \le \dots\}$35, with $\{T_1 \le T_2 \le \dots\}$36 denoting an active bursty state and $\{T_1 \le T_2 \le \dots\}$37 an inactive sparse state. The intensity becomes

$\{T_1 \le T_2 \le \dots\}$38

with

$\{T_1 \le T_2 \le \dots\}$39

and rank-based transition rates

$\{T_1 \le T_2 \le \dots\}$40

The stationary active-state probability is

$\{T_1 \le T_2 \le \dots\}$41

so dyads with large $\{T_1 \le T_2 \le \dots\}$42 spend more time in the active state (Ward et al., 2020).

This rank-structured hierarchy encodes several behavioral phenomena simultaneously. The Hawkes term models the winner effect, because recent wins elevate future win intensity. The active state produces bursting, since while $\{T_1 \le T_2 \le \dots\}$43 the dyad follows Hawkes dynamics rather than a low-rate Poisson background. Pair-flips are permitted because subordinate-to-dominant excitation is small but nonzero, especially when ranks are close, and because the latent state process allows noisy inactive periods (Ward et al., 2020).

On simulated data, C-MMHP and C-DCHP recovered the true latent ranks well, while C-HP showed noticeable bias for some nodes. On real mice-aggression data, C-MMHP, C-DCHP, and C-HP gave broadly consistent rank estimates, but C-MMHP often had narrower credible intervals and slightly best overall predictive performance; active-state interactions exhibited a more linear hierarchy, whereas inactive-state interactions were noisier (Ward et al., 2020).

Beyond social interaction networks, hierarchical Hawkes structure is natural in seismology, finance, crime, and epidemics. ETAS already combines background events, magnitude-dependent offspring productivity, and spatial kernels; multivariate financial Hawkes models organize excitation across assets, sectors, or markets; crime and epidemic models introduce spatial and demographic nesting (Laub et al., 2024, Lima, 2020). The central limitations are likewise recurrent: identifiability of background versus excitation, parameter growth of order $\{T_1 \le T_2 \le \dots\}$44 in large multivariate systems, slow posterior computation in hierarchical Bayesian formulations, partial or binned observations, and, in segmentation-based approaches, deterministic boundaries with no quantified uncertainty (Laub et al., 2024, Zhou et al., 2019). These constraints explain why hierarchical Hawkes research has focused on shrinkage, low-rank structure, structured regularization, cumulant-based estimation, and latent-variable augmentation rather than on unrestricted parameterization.