Skip to content
All library documents

Attempting Multivariate Hawkes Simulation with Ogata Thinning

Article Quant Q&A · Author: jraffaud

Summary

The document presents a Python attempt to simulate a multivariate Hawkes process using Ogata’s thinning algorithm. It represents event times for each component as separate arrays and computes component intensities from baseline rates plus exponentially decaying contributions from past events. The proposed loop samples a candidate waiting time from the summed intensity, selects a component, and appends the event time to that component’s history.

The questioner reports that the maximum intensity grows unexpectedly and that the simulation does not match a reference. The document provides no answer diagnosing the issue or validating the implementation. It is useful as a setup for understanding a simulation problem, but readers should not treat the code as a verified algorithm. In particular, the text gives no results, checks, or discussion of assumptions needed to assess correctness or stability.

Key ideas

  • The code represents a multivariate process with a separate event-time history for each component.
  • Each component’s intensity combines a baseline with exponential contributions from earlier events.
  • Ogata thinning is used to propose event times based on the summed intensity.
  • The reported implementation produces unexpectedly large intensities and is not validated in the document.

Tags

Full text
# Multivariate Hawkes Process Simulation


# Multivariate Hawkes Process Simulation












I am trying to implement Ogata's thinning algorithm to simulate multivariate Hawkes Processes in Python (the algorithm can be found here: https://www.math.fsu.edu/~ychen/research/Thinning%20algorithm.pdf), but I'm running into some trouble.

Here's my full code so far:

```
def multivariate_cif(t, T, mu, alpha, beta):

    n = len(mu)

    for i in range(n):
        for j in range(n):
            mu[i] += sum(alpha[i, j] * 
                         np.exp(-beta[i, j] * 
                         np.subtract(t, T[j][np.where(T[j]<t)])))

    return mu

def multivariate_simulation(time, mu, alpha, beta):

    T = [np.array([]) for _ in range(len(mu))]
    n = np.zeros((len(mu)))
    b = 0; s = 0

    while (b < time):

        M = sum(multivariate_cif(s, T, mu, alpha, beta))
        s += -np.log(np.random.uniform(0,1))/M

        D = np.random.uniform(0, M)

        if (D <= sum(multivariate_cif(s, T, mu, alpha, beta))):

            k = 0

            while (D > sum(multivariate_cif(s, T, mu, alpha, beta)[:k+1])):

                k += 1

            n[k] += 1
            T[k] = np.append(T[k], s)

        b +=1

    return T
```

Here T is a list of arrays, with each array containing arrival times for one of the counting processes.

Basically when I run this simulation my max value "M" explodes, I can't replicate the reference simulation for the life of me despite having each line down exactly like the ref.

Anyone who has successfully run Ogata simulation, could you please shed some light on where I'm getting lost?

Thanks a lot!!

Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)

This summary was written by Stratmill's research agent from the original; it is not a copy of the source.