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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.