Calculating a correlation integral with np.random

I'm trying to calculate a correlation integral using Python, to estimate a parameter which I need to solve a PDE via finite element. The correlation integral reads

where omega and c-bar are positive constants. The m terms are N(0,1) Gaussian random variables. The z parameter is the length coordinate. In discrete form it can be written

Here lowercase m is a 1xM vector containing N(0,1) Gaussian random variables. I want to model the case when the z values are drawn according to a Poisson distribution. I achieve this by using a cumulative sum.

The theory says that this correlation integral should be strictly positive as it is related to the power spectral density. However, when I run this a bunch of times, the positive and negative results for the correlation integral are 50-50. Can anyone shed some light on why it's a 50-50 split? My code (here the upper limits on both sums N=M):

import matplotlib.pyplot as plt
import numpy as np
import random
#variance in N(0,1) rvs
sigma = 1
#Parameters
omega = 5*10**6
lambda3 = 0.005
cbar = omega*lambda3/(2*np.pi)
#Length of vectors
M = 100
#Create Poisson coordinate vector
lam = 4
poisson_z = np.random.poisson(lam,M)
#Poisson z coordinate vector
z_vec_poisson = np.cumsum(poisson_z)
#Poisson increment vector
poisson_diff_vec = abs(np.insert(np.diff(poisson_z),0,0))

#Empty list to store correlation integrals
corr_full_cos = []
#Calculate correlation integral 'runs' times
runs = 10000
for k in range(1,runs):
    #Percentage complete
    print(round(k/(runs-1),3)*100,'%')

    #create empty vector for expectation sum
    summy=np.zeros(M)
    for p in range(1,M):
        #Draw p N(0,1) random vectors
        m_vector_p  = np.random.normal(0,sigma,M)
        #Update for each sum +=
        summy += m_vector_p[0]*m_vector_p
        #Hadamard multiply summy with poisson increment vector
        summy_delta_poisson_z = (1/M)*np.multiply(summy,poisson_diff_vec)
    
    
    #Evaluate correlation integral for Poisson z
    corr_int_general_cos = 2*np.dot(summy_delta_poisson_z,np.cos(2*omega*z_vec_poisson/cbar))

    #Store correlation integral value 
    corr_full_cos.append(corr_int_general_cos)

#Negative and positive correlation integrals out of k runs
pos_count,neg_count = 0,0
for num in  corr_full_cos: 
      
    # checking condition 
    if num >= 0: 
        pos_count += 1
  
    else: 
        neg_count += 1
print("Positive : ", (pos_count/runs)*100,'%') 
print("Negative : ", (neg_count/runs)*100,'%')