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

[Francesco Potorti`]

# Put in the public domain by:
# Francesco Potorti` <address@hidden> 1995

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