def divisorGenerator(n):
    large_divisors = []
    for i in range(1, int(math.sqrt(n) + 1)):
        if n % i == 0:
            yield i
            if i*i != n:
                large_divisors.append(n / i)
    for divisor in reversed(large_divisors):
        yield divisor

def change_dataLength(vector):
    '''
    Function that removes last observations from input series, such that the data
    can be split in more than 2 ranges.
    '''
    found = False
    n= len(vector)
    while found==False:
        full_divisors = list(divisorGenerator(n))
        if len(full_divisors) <= 3:
            n -= 1
        else:
            found=True
    new_vector = vector[0:n]
    full_divisors=full_divisors[:-1] #last divisor is the length of the series, not useful for ranges
    return new_vector, full_divisors

def RS_function(input_series):
    incs = input_series
    mean_inc = np.sum(incs) / len(incs)
    deviations = incs - mean_inc
    Z = np.cumsum(deviations)
    R = max(Z) - min(Z)
    S = np.std(incs, ddof=1)
    return R/S

def calc_Hurst(original_series):
    input_series, splits = change_dataLength(original_series)
    window_sizes = len(input_series)/np.asarray(splits)
    window_sizes = window_sizes.astype(int)
    if np.sum(window_sizes < 4)>0:
         window_sizes = window_sizes[window_sizes >= 4]


    err = np.geterr()
    np.seterr(all='raise')

    R_S = []
    for w in window_sizes:
        rs = []
        for start in range(0, len(input_series), w):
            if (start+w)>len(input_series):
                break
            temp_rs = RS_function(input_series[start:start+w])
            if temp_rs != 0:
                rs.append(temp_rs)
        R_S.append(np.mean(rs))

    A = np.vstack([np.log10(window_sizes), np.ones(len(R_S))]).T
    H, c = np.linalg.lstsq(A, np.log10(R_S), rcond=-1)[0]
    np.seterr(**err)

    c = 10**c
    return H, c, [window_sizes, R_S]

start= 0
finish = 144000
delta = 0.25


times = np.arange(start,finish, delta)
B = np.zeros(len(times))
for i in range(len(times)):
    if i == 0:
        B[i] = 0
else:
    u = np.random.normal(0,1)
    v = u * np.sqrt(times[i] - times[i-1])
    B[i]= B[i-1]+ v