function [H, r, x, y] = hurst3(A)
% [H,r,x,y] = hurst3(A)
% Computes the Hurst parameter H of the vector A using the
% variance-time plot method.  In this version A is a walk,
% i.e. the cumulative sum of a noise a.  Modified from an
% Octave program with original comments:
% Uses the variance computed over intervals of different lengths.
% Plotted on a log-log diagram, var (x) versus the interval length
% (y) should be least-squares interpolated by a line of slope
% beta=2H.  r is the correlation coefficient and gives a
% "reliability factor" of the H estimate: values higher than
% 0.9 should be expected.
%
% I don't know about any data on the precision of this method.
% All I know is that it cannot provide a confidence interval measure.
%
% Leland, Taqqu, Willinger, Wilson, "On the Self-Similar Nature
% of Ethernet Traffic (Extended Version)", IEEE/ACM Trans. on
% Networking Vol.2, Num.1, 1994.

% Put in the public domain by Francesco Potorti`
% <pot@cnuce.cnr.it> Mon Aug 28 17:07:26 MET 1995.
% Modified for MATLAB by Gerry Middleton, December 1995.
% Uses function tf = isvec(A) to test that A is a vector.
% This version reduces the minimum no of data required, and
% changes the input from a noise to a walk. In general, the
% more data the better but if you have less than 631 points
% use this version.     It still requires > 398 points.

% ts = isvec(A);  % test that A is a vector

% m is the minimum interval length over which the var 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.
% maxp is the maximum number of points that must be used.

m = 10;
M = 10;
k = 5;
minp = 4;        % originally 5, used only to calculate minr
maxp = 20;

r = length(A);    % test for enough data
minr = ceil(m*M*10^((minp-1)/k));
if (r < minr)
  error (sprintf ('A should have at least %d elements\n', minr));
end

% n is the number of points we use to interpolate the line on a
% log-log plot whose slope is 2H.  x is a vector containing
% their abscissae, y contains their ordinates, that is, the
% variance of the increments of the vector A at distances x.

n = min(maxp, floor(k * log10(r/m/M)));
x = floor(logspace (log10(m), log10(r/M), n))';

% now calculate y

cumA = A;        % if data are noise, change to cumA = cumsum(A);
for i = 1:n
  xi = x(i);
  Y = cumA(xi:xi:r) - [0; cumA(xi:xi:r-xi)];
  sets = length(Y);
  y(i) = (sum(Y.*Y) - cumA(sets*xi)^2/sets) / (sets-1);
end

% take logs and convert y to a column vector

x = log(x);
y = log(y');

% solve for H -- alfabeta(2) is the slope of the log plot
I = ones(size(x));
alfabeta = [I, x] \ y;
H = alfabeta(2) / 2;
if (nargout > 1)
   r = corrcoef(x, y);
end

% plot

plot(x,y, x,y,'o');
title('log(variance at delta x) vs log(delta x)');
xlabel('log(delta x)'), ylabel('log(variance)');