I write a function which takes as input 2 arrays of zeros and ones ~8000 elements per array. Input array size is variable, can have much higher number of elements. The density of ones is expected to be higher than zeros for most cases, but outliers may exist. My function `eps` calculates a statistic on these array pairs and returns the output. This statistic is just a way to find if 2 timeseries data correlate with each other. 
 Relatively trivial operations of checking for 0 and noting the index where 0 is found in array and then some calculations. I tried my best to optimize for speed but the best I could get is 4.5 ~5 seconds (for 18k array pairs) using `timeit` (1000 runs) library. Fast execution is important as I need to run this function on billions of array pairs.

UPDATE: The solution I tried from the answers below improved runtime for larger inputs than my original approach. However, I still face the problem of running this statistical function a 1000 times per array pair. The background for this is to perform a permutation test on the statistic with 'x' (1000) number of permutations to arrive at a significance level at p= .01
Some asked below why I need to run this function billions of times. So my data is a timeseries data for events: 'E1,E2,E3.....En' I need to find combinations of all possible event pairs from this data i.e. nCr. (E1,E2), (E2,E3)....nCr times. This leads to a large number of event pairs about billions of event pairs for my data. Add to that the time complexity of a permutation test (1000x times/pair). Any suggestions if any other parts like permutation tests can be improved. Thanks  
UPDATE2: Is size always ~8000? No it depends on input data, array size can be between 2000 ~ 10,000.
What's your tau value? tau is a random positive integer or float and we jus test for different tau values, so I cannot given a certain Tau.

    #e.g. inputs
    #ts_1 = [0,1,1,0,0,1,1,0,......] #timeseries1
    #ts_2 = [1,1,1,1,1,1,1,0,......] #timeseries2
    # tau = positive integer/float   #time lag

    import numpy as np
    from itertools import product

     def eps(ts_1, ts_2, tau):  
        Q_tau = 0
        q_tau = 0
        
        event_index1, = np.where(np.array(ts_1) == 0)
        n1 = event_index1.shape[0]
        
        event_index2, = np.where(np.array(ts_2) == 0)
        n2 = event_index2.shape[0]
        
        if (n1 != 0 and n2 != 0):
            matching_idx = set(event_index1).intersection(event_index2)
            c_ij = c_ji = 0.5 *len(matching_idx)
            
            for x,y in product(event_index1,event_index2):
                if x-y > 0 and (x-y)<= tau:
                    c_ij += 1
                elif y-x > 0 and (y-x) <= tau:
                    c_ji += 1                     
                      
            Q_tau = (c_ij+c_ji)/math.sqrt( n1 * n2 )
            q_tau = (c_ij - c_ji)/math.sqrt( n1 * n2 )
        
        return Q_tau, q_tau

    #permutation test
    
     def eps_permutationtest(ts_1,ts_2,permutations=1000,Q_tau): #Q_tau= original statistic to evaluate significance for
            
         lst_Q_tau = []    # list to hold Q_tau values from each shuffle 
            
         for i in range(0,permutations):
             np.random.shuffle(ts_1)         #shuffle ts_1 everytime
             Q_tau = eps(ts_1 = ts_1, ts_2 = ts_2 , tau= 5) 
             lst_Q_tau.append(Q_tau)
         lst_Q_tau= np.array(lst_Q_tau)
            
         p_val = len(np.where(lst_Q_tau>= [Q_tau])[0])/permutations
            
            
         return p_val  # significant if p < .01 alpha