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`
% 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)');