This is an implementation of the blocking method of Jonsson (PHYSICAL REVIEW E 98, 043304 (2018)) in Julia. It takes an array as input and will output the estimated error of the mean value using rigorous blocking arguments.

blockjonsson.jl

function blockjonsson(a::AbstractArray)
    d = log2(length(a))
    x = copy(a)
    if (d - floor(d)) != 0
        @info "Warning: Data size = $(floor(2^d)), is not a power of 2."
        @info "Truncating data to $(2^floor(d))."
        x = x[1:(2^round(Int, floor(d)))]
    end
    d = round(Int, floor(d))
    n = 2^d
    s = zeros(d)
    gamma = zeros(d)
    mu = mean(x)

    for i in 1:d
        n = length(x)
        gamma[i] = sum((x[1:(n - 1)] .- mu) .* (x[2:n] .- mu)) / n
        s[i] = var(x)
        x = (x[1:2:end] .+ x[2:2:end]) ./ 2.0
    end

    # Observator
    range_d = reverse!(collect(1:d))
    twopowers = fill(2.0, length(range_d))
    M = reverse!(cumsum(reverse!((gamma ./ s) .^ 2 .* (twopowers .^ range_d))))

    k_val = 0
    for k in 1:d
        if M[k] < q[k]
            k_val = k - 1
            break
        end
    end

    if k_val >= d - 1
        println("Warning: Use more data")
    end

    sem = sqrt(s[k_val + 1] / 2^(d - k_val))

    return sem
end