📊 statistics

Parametric inference for the discretely observed multivariate Hawkes process using particle Markov Chain Monte Carlo

This paper proposes a particle Markov Chain Monte Carlo method that utilizes sequential Monte Carlo to estimate the intractable likelihood of discretely observed multivariate Hawkes processes, demonstrating superior performance over existing approaches on both simulated data and real-world terrorist activity records.

Original authors: Jason J. Lambe, Feng Chen, Tom Stindl, Tsz-Kit Jeffrey Kwan

Published 2026-06-18

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Jason J. Lambe, Feng Chen, Tom Stindl, Tsz-Kit Jeffrey Kwan

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to understand how a rumor spreads through a school. You know that when one student hears a rumor, they are likely to tell a few friends, who then tell more friends. This is a "self-exciting" process: one event triggers more events.

Now, imagine you are the principal, but you don't get a live feed of every single conversation. Instead, you only get a daily report at the end of the day saying, "Today, 5 rumors started in the cafeteria and 3 in the library." You don't know when exactly they happened or who told whom.

This is the problem the paper solves. It deals with a complex mathematical model called a Multivariate Hawkes Process (which tracks multiple types of events, like rumors in different locations) when the data is "discretely observed" (only daily counts, not exact times).

Here is a breakdown of what the authors did, using simple analogies:

1. The Problem: The "Blind" Likelihood

In statistics, to figure out the rules of how the rumors spread (the parameters), you usually need to calculate something called a "likelihood." Think of the likelihood as a scorecard that tells you how well your theory fits the data.

The Catch: When you only have daily totals (like "5 rumors today") instead of exact times, the math to calculate this scorecard becomes impossible to solve directly. It's like trying to solve a puzzle where half the pieces are missing and the picture on the box is blurry. Traditional methods fail here.

2. The Solution: The "Guess-and-Check" Team (SMC)

The authors created a clever way to estimate that impossible scorecard. They used a method called Sequential Monte Carlo (SMC).

The Analogy: Imagine you have a team of 1,000 detectives (called "particles"). Each detective tries to guess what the hidden conversations looked like to result in the daily count of 5 rumors.
- Some detectives guess the rumors happened at 9:00 AM, others at 11:00 AM.
- They check their guesses against the rules of the school (the math model).
- If a detective's guess makes no sense (e.g., they guessed 10 rumors happened when the report said 5), they are eliminated.
- The detectives who have good guesses get "more votes" (weights).
The Innovation: The authors realized that if they just let the detectives guess randomly, many would waste time on impossible scenarios. So, they designed a specific rule for the detectives: "If the report says 5 rumors happened, your 5 guesses must be spread out evenly across the day." This simple trick (called an "ordered uniform proposal") made the team of detectives much smarter and faster, ensuring they didn't waste time on bad guesses.

3. The Engine: The "Fair Coin" Machine (PMMH)

Once the team of detectives gives them an estimated scorecard, the authors need to find the best set of rules for the school. They use a method called Pseudo-Marginal Metropolis-Hastings (PMMH).

The Analogy: Imagine you are trying to find the highest peak in a foggy mountain range (the best parameters). You can't see the whole mountain.
- You take a step to a new spot.
- You ask your team of detectives to estimate the height of the new spot.
- Because the team is unbiased (they don't systematically lie about the height), you can trust their estimate to decide whether to stay at the new spot or go back.
- Over time, by taking many steps and listening to the team, you map out the entire mountain range and find the highest peak.
Why it matters: This allows them to not just find the "best" answer, but also to know how confident they are in that answer (like saying, "We are 95% sure the peak is between these two points").

4. The Results: Better than the Competition

The authors tested their method against other existing ways of solving this problem using simulated data (fake school rumor scenarios).

Accuracy: Their method was more accurate (lower error) than the other methods.
Speed: It was significantly faster. While other methods might need a huge team of detectives or take hours to run, their method got great results with a smaller team and less time.
Flexibility: It works even if the daily reports come in at different times (e.g., some days are 12 hours, some are 24 hours) or if the background noise of the school changes throughout the day.

5. Real-World Test: Terrorist Activity in Afghanistan and Pakistan

To prove it works in the real world, they applied their method to real data: daily counts of terrorist attacks in two neighboring regions (Kabul/Nangarhar in Afghanistan and Khyber Pakhtunkhwa in Pakistan) from 2018 to 2021.

What they found:
- Attacks in one region tend to trigger attacks in the same region (self-excitation).
- Crucially: Attacks in Pakistan (Khyber Pakhtunkhwa) strongly triggered attacks in Afghanistan (Kabul/Nangarhar).
- They found that while attacks within a single region happen very quickly (within hours), the "cross-border" trigger takes about a month to play out.
- They also noticed a spike in attacks in Afghanistan later in the timeline, which aligned with the period of the US military withdrawal.

Summary

The paper introduces a new statistical toolkit that allows researchers to understand complex, "contagious" events (like rumors, financial trades, or terrorist attacks) even when they only have rough, daily summaries of the data. By using a smart team of "detectives" to guess the missing details and a "fair coin" machine to find the best rules, they can get accurate answers faster and more reliably than previous methods.

Technical Summary: Parametric Inference for Discretely Observed Multivariate Hawkes Processes

Problem Statement
The Multivariate Hawkes Process (MHP) is a standard model for analyzing multidimensional event sequences characterized by self-excitation and cross-excitation. While Maximum Likelihood Estimation (MLE) and Expectation Maximization (EM) algorithms are well-established for MHPs with precise event times, they fail when data is observed only in discrete intervals (aggregated counts). In such scenarios, the likelihood function becomes analytically intractable because the exact event times and types within the intervals are unobserved latent variables. Existing approaches to this problem, such as Integer-Valued Autoregressive (INAR) approximations, modified Monte Carlo EM (MCEM), coarse-grained Hawkes processes, and Bayesian MCMC, often suffer from bias, computational inefficiency, or an inability to handle non-stationary background rates and unequal censoring intervals.

Methodology
The authors propose a framework for parametric inference using Pseudo-Marginal Metropolis-Hastings (PMMH) combined with Sequential Monte Carlo (SMC).

State-Space Representation: The unobserved event times and types within each observation interval are treated as latent variables in a state-space model. The likelihood of the aggregated data is factorized into a product of conditional probabilities over time intervals.
Unbiased Likelihood Estimation via SMC: To address the intractable likelihood, the authors employ an SMC algorithm to generate an unbiased estimate of the likelihood function.
- Proposal Distribution: A key innovation is the use of an ordered uniform proposal distribution for latent event times within an interval, conditioned on the observed counts. This ensures that every proposed particle is consistent with the observed data (unlike Poisson proposals which may generate particles with zero weight). Event types are assigned via a multivariate hypergeometric distribution.
- Adaptive Resampling: The SMC algorithm utilizes adaptive resampling based on the Effective Sample Size (ESS) to mitigate particle degeneracy while minimizing the introduction of Monte Carlo noise.
PMMH Algorithm: The unbiased SMC likelihood estimate is substituted for the true likelihood within a Metropolis-Hastings algorithm. Theoretical results (Andrieu and Roberts, 2009) guarantee that the stationary distribution of the resulting Markov chain is the true posterior distribution of the parameters, provided the likelihood estimator is unbiased.
Parameter Estimation: The posterior mean or median serves as the point estimate, with credible intervals derived from the MCMC sample quantiles. The method supports time-varying background rates (modeled via B-splines) and handles unequally sized observation intervals.

Key Contributions

Unbiased Inference: The paper provides a method to perform exact Bayesian inference (in the limit of infinite particles) for discretely observed MHPs by constructing an unbiased likelihood estimator.
Efficient Proposal Mechanism: The introduction of the ordered uniform proposal significantly reduces the variance of the likelihood estimate compared to previous Poisson-based proposals. This allows the PMMH algorithm to converge with a much smaller number of particles (orders of magnitude fewer), drastically improving computational efficiency.
Flexibility: The framework accommodates non-stationary background rates and irregular observation windows without requiring data augmentation or model discretization.

Results

Simulation Studies:
- Validity: SMC estimates of event probabilities converge to true Monte Carlo estimates with low bias and variance as the number of particles increases.
- Proposal Efficiency: The ordered uniform proposal yields significantly lower variance in log-likelihood estimates compared to Poisson proposals, enabling the use of fewer particles (e.g., $J=10$ ) while maintaining accuracy.
- Performance Comparison: In simulations across varying aggregation levels ( $\Delta \in \{0.1, 0.5, 1.0\}$ ), the PMMH estimator demonstrated minimal bias and well-calibrated confidence intervals. It outperformed MCEM and Coarse-Grained Hawkes Process (CGHP) methods in terms of Mean Squared Error (MSE), particularly at higher aggregation levels where those methods exhibited significant bias.
- Computational Speed: The PMMH method was faster than the Bayesian MCMC approach (which could not be parallelized in the referenced implementation) and the MCEM/CGHP methods, scaling efficiently with parallel computation of particles.
Empirical Application:
- The method was applied to daily terror attack counts in Kabul–Nangarhar (Afghanistan) and Khyber Pakhtunkhwa (Pakistan) from 2018 to 2021.
- The model successfully captured the dynamics, revealing significant self-excitation in both regions and strong cross-excitation from Pakistan to Afghanistan ( $\hat{\eta}_{1,2} = 0.438$ ).
- The estimated background rates reflected the acceleration of attacks in Afghanistan following the US withdrawal, consistent with external reports. The model fit the observed data well, with empirical counts falling within the simulated 95% quantiles.

Significance and Claims
The authors claim that their method offers a statistically rigorous and computationally efficient alternative for analyzing discretely observed MHPs. By leveraging the unbiased nature of SMC within a PMMH framework, the method avoids the bias inherent in MCEM and the approximation errors of coarse-grained models. The specific contribution of the ordered uniform proposal is highlighted as the critical factor enabling the method's speed and stability, making it competitive with high-quality alternatives while providing immediate uncertainty quantification (credible intervals) and flexibility for non-stationary processes. The paper modestly notes that while the estimators perform well empirically, the theoretical proof of consistency and asymptotic normality for the specific MHP context remains a subject for future work, relying currently on general SMC/PMMH theory.