for x,y in product(event_index1,event_index2) looks like an efficiency killer. If the lengths of event_index1, event_index2 are L1, L2 respectively, this loop has O(L1 * L2) time complexity. You may get away with O(L1 + log L2) (or O(L2 + log L1), depending on which one is larger.

Notice that for any given x the c_ij is incremented by the number of ys in the [x - tau, x)range. Two binary searches over event_index2 would give you this range, along the lines of

    for x in event_index1:
       c_ij += bisect_left(event_index2, x) - bisect_left(event_index2, x - tau)

and ditto for c_ji. Of course, from bisect import bisect_left.

for x,y in product(event_index1,event_index2) looks like an efficiency killer. If the lengths of event_index1, event_index2 are L1, L2 respectively, this loop has O(L1 * L2) time complexity. You may get away with O(L1 * log L2) (or O(L2 * log L1), depending on which one is larger.

Notice that for any given x the c_ij is incremented by the number of ys in the [x - tau, x)range. Two binary searches over event_index2 would give you this range, along the lines of

    for x in event_index1:
       c_ij += bisect_left(event_index2, x) - bisect_left(event_index2, x - tau)

and ditto for c_ji. Of course, from bisect import bisect_left.