help-octave
[Top][All Lists]

## hurst.m: compute the Hurst parameter of an array of samples

 From: Francesco Potorti` Subject: hurst.m: compute the Hurst parameter of an array of samples Date: Thu, 20 Jul 1995 13:20 +0100 (MET)

# Put in the public domain by:

function [H, r, IDC, len] = hurst (A)

# Computes the Hurst parameter H of the vector column A.
#
# Use the Index of Dispersion for Counts (IDC, or peakedness), that is, the
# ratio of the variance over mean, computed over intervals of different
# lenghts.  Plotted on a log-log diagram, IDC versus the interval length len
# should be least-squares interpolated by a line of slope beta=2H-1.  r is
# the correlation coefficient and gives a "reliability factor" of the H
# estimate: values higher than 0.9 should be expected.
#
# Leland, Taqqu, Willinger, Wilson, "On the Self-Similar Nature of Ethernet
# Traffic (Extended Version)", IEEE/ACM Trans. on Networking Vol.2, Num.1,
# 1994.

if (!is_vector(A) || columns(A) != 1)
error ("A should be a column vector\n");
endif

# m is the minimum interval length over which the IDC is computed.  M is
# the minimum number of intervals over which that same quantity is
# computed.  k is the number of interval lengths per decade used for the
# computation.  minp is the minimum number of points that must be used.
m = 3;
M = 10;
k = 8;
minp = 5;

r = rows (A);
minr = ceil (m*M*10^((minp-1)/k));
if (r < minr)
error (sprintf ("The argument should have at least %d rows\n", minr));
endif
n = floor (k * log10(r/m/M));
IDC = floor (logspace (log10(m), log10(r/M), n))';
for i = 1:n
sets = floor (r/IDC(i));
Len = sum (reshape (postpad (A, sets*IDC(i)), IDC(i), sets));
s = sum (Len);
len(i) = (sumsq(Len)*sets/s - s) / (sets-1);
endfor
logIDC = log (IDC);
loglen = log (len);
alfabeta = [ones (size (IDC)), logIDC] \ loglen;
H = (alfabeta(2)+1)/2;
if (nargout > 1) r = corrcoef (logIDC, loglen); endif

endfunction