def f(z):
    return ((1-z)*np.exp(-d/z))/(((1-z)**(1+d))*(z**(1-d)))
lhs=integrate.quad(f,0,0.5)
def rhs(p):
    return integrate.quad(-f,0.5, p)
p_star= fsolve(rhs-lhs,0.75) # will depend on time series only indrectly when d will be optimized
for i in np.arange(N):
    z=Z[i,:]
    main = np.exp (kz+m*np.arange(w))
    cumsum_t=np.cumsum(1/main)
    final_sum= main*cumsum_t
    t_solution= # index i where final_sum[i]> p_star/(1-p_star)) # not implemented yet