Matlab:尝试通过直方图盒计数从时间序列估计多重分形谱

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

我正在使用耶鲁分形页面中的方法: http://classes.yale.edu/fractals/MultiFractals/Moments/TSMoments/TSMoments.html
这组讲座幻灯片(幻灯片 32)也对此进行了阐述: http://multiscale.emsl.pnl.gov/docs/multifractal.pdf

这个想法是你得到一个数据集,并通过许多直方图检查它,增加条的数量,即分辨率。一旦分辨率足够高,一些条形就会取零值。在每个阶段,我们取落入每个直方图箱(忽略任何零值箱)的结果数量,将其除以数据集的总大小,并将其提高到幂 q。此操作给出了给定时刻和 bin 大小的 'partition function'。引用上面链接的教程:“提供了对非均匀性的选择性表征 测量,正 q 突出最密集的区域,负 q 突出最平滑的区域。

所以我在 Matlab 中使用直方图函数,遍历 bin 大小,对所有非零 bin 内容求和,等等。但是我的分区函数输出数组只是一堆 1。我看不出出了什么问题,其他人可以吗?

英特尔、思科、苹果和其他公司的数据可在同一耶鲁网站上获得:yale.edu/fractals/MultiFractals/Finf(a)/Finf(a).html

N.B。 intel指的是我原先用作数据集的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

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

您的代码似乎总体上没有错误,但我做了一些更改,因为您在循环中执行了不必要的重复(我将外循环移到了内部,"vectorized" 因为所有时刻的计算都可以同时执行给定直方图。此外,它正在构建花费时间最长的直方图)。

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) ....