# 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!!