Matlab: trying to estimate multifractal spectrum from time series by histogram box-counting

Question

I am using the approach from this Yale page on fractals: http://classes.yale.edu/fractals/MultiFractals/Moments/TSMoments/TSMoments.html
which is also expounded on this set of lecture slides (slide 32): http://multiscale.emsl.pnl.gov/docs/multifractal.pdf

The idea is that you get a dataset, and examine it through many histograms with increasing numbers of bars i.e. resolution. Once resolution is high enough, some bars take the value zero. At each of these stages, we take the number of results that fall into each histogram bin (neglecting any zero-valued bins), divide it by the total size of the dataset, and raise this to a power, q. This operation gives the 'partition function' for a given moment and bin size. Quoting the above linked tutorial: "Provides a selective characterization of the nonhomogeneity of the measure, positive q’s accentuating the densest regions and negative q’s the smoothest regions."

So I'm using the histogram function in Matlab, looping over bin sizes, summing over all the non-zero bin contents, and so forth. But my output array of partition functions is just a bunch of 1s. I can't see what's going wrong, can anybody else?

Data for intel, cisco, apple and others is available on the same Yale website: yale.edu/fractals/MultiFractals/Finf(a)/Finf(a).html

N.B. intel refers to the intel stock price I was originally using as the dataset.

lower = 1; %set lowest level of histogram resolution (bin size)
upper = 300; %set highest level of histogram resolution (bin size)

qlow = -20; %set lowest moment exponent
qhigh = 20; %set highet moment exponent
qstep = 0.25; %set step size between moment exponents

qn= ((qhigh-qlow)/qstep) + 1; %calculates number of steps given qlow, qhigh, qstep

qvalues= linspace(qlow, qhigh, qn); %creates a vector of q values given above parameters

m = min(intel); %find the maximum of the dataset
M = max(intel); %find the minimum of the dataset

for Q = 1:length(qvalues) %loop over moment exponents q 


    for k = lower:upper %loop over bin sizes

        counts = hist(intel, k); %unpack all k histogram height values into 'counts'
        counts(counts==0) = []; %delete all zero values in ''counts

        Zq = counts ./ length(intel);
        Zq = Zq .^ qvalues(Q);
        Zq = sum(Zq);

        partitions(k-(lower-1), Q) = Zq; %store Zq in the kth row and the Qth column of 'partitions'

    end

end

Answer 1

Your code seems to be generally bug-free but I made some changes since you perform needless repetitions over loops (I moved the outer loop inside and "vectorized" it since all moment calculations can be performed simultaneously for a given histogram. Also, it is building the histogram that takes longest).

intel = cumsum(randn(64,1)); % <-- mock random walk
Ni =length(intel);

%figure, plot(intel)


lower = 1; %set lowest level of histogram resolution (bin size)
upper = 300; %set highest level of histogram resolution (bin size)

qlow = -20; %set lowest moment exponent
qhigh = 20; %set highet moment exponent
qstep = 0.25; %set step size between moment exponents

qn= ((qhigh-qlow)/qstep) + 1; %calculates number of steps given qlow, qhigh, qstep

qvalues= linspace(qlow, qhigh, qn); %creates a vector of q values given above parameters

m = min(intel); %find the maximum of the dataset
M = max(intel); %find the minimum of the dataset

partitions = zeros(upper-lower+1,length(qvalues));
for k = lower:upper %loop over bin sizes

    % (1) Select a bin size r and partition [m,M] into intervals of size r:
        % [m, m+r), [m+r, m+2r), ..., [m+kr, M], where m+kr < M <= m+(k+1)r.
    % Call these bins B0, ..., Bk.

    edges = linspace(m,M,k+1);
    edges(end)=Inf;

    % (2) For each j, 0 <= j <= k, count the number of xi that lie in bin Bj. Call this number nj. Ignore all nj that equal 0 after all the xi have been counted..

    counts = histc(intel, edges); %unpack all k histogram height values into 'counts'
    counts(counts==0) = []; %delete all zero values in ''counts

    % (3) Now compute the qth moment, Mrq = (n0/N)q + ... + (nk/N)q, where the sum is over all nonzero ni.

    % Zq = counts/Ni;

    partitions(k, :) = sum( (counts/Ni) .^ qvalues); %store Zq in the kth row and the Qth column of 'partitions'

end

figure, hold on
loglog(1./[1:k]', partitions(:,1),'g.-')
loglog(1./[1:k]', partitions(:,80),'b.-')
loglog(1./[1:k]', partitions(:,160),'r.-')

% (4) Perform linear regressions here to get alpha(r) ....