% The Hurst exponent

%--------------------------------------------------------------------------

% The first 20 lines of code are a small test driver.

% You can delete or comment out this part when you are done validating the

% function to your satisfaction.

%

% Bill Davidson, quellen@yahoo.com

% 13 Nov 2005

function []=hurst_exponent()

disp('testing Hurst calculation');

n=100;

data=rand(1,n);

plot(data);

hurst=estimate_hurst_exponent(data);

[s,err]=sprintf('Hurst exponent = %.2f',hurst);disp(s);

%--------------------------------------------------------------------------

% This function does dispersional analysis on a data series, then does a

% Matlab polyfit to a log-log plot to estimate the Hurst exponent of the

% series.

%

% This algorithm is far faster than a full-blown implementation of Hurst's

% algorithm.

I got the idea from a 2000 PhD dissertation by Hendrik J

% Blok, and I make no guarantees whatsoever about the rigor of this approach

% or the accuracy of results.

Use it at your own risk.

%

% Bill Davidson

% 21 Oct 2003

function [hurst] = estimate_hurst_exponent(data0)

% data set

data=data0;  

% make a local copy

[M,npoints]=size(data0);

yvals=zeros(1,npoints);

xvals=zeros(1,npoints);

data2=zeros(1,npoints);

index=0;

binsize=1;

while npoints>4

y=std(data);

index=index+1;

xvals(index)=binsize;

yvals(index)=binsize*y;

npoints=fix(npoints/2);

binsize=binsize*2;

for ipoints=1:npoints % average adjacent points in pairs

data2(ipoints)=(data(2*ipoints)+data((2*ipoints)-1))*0.5;

end

data=data2(1:npoints);

end % while

xvals=xvals(1:index);

yvals=yvals(1:index);

logx=log(xvals);

logy=log(yvals);

p2=polyfit(logx,logy,1);

hurst=p2(1); % Hurst exponent is the slope of the linear fit of log-log plot

return;

function get_hits(fname, levels, delta, deleteFirst, isTimeSeries, isIncremental, y_eps)

% hitting times, points and subcrossings for continuous process
%
% get_hits(fname, levels, delta, deleteFirst, isTimeSeries, isIncremental, y_eps)
%
% fname: input file name without extension, x(t(n)) and t(n) in files fname.dat and fname.tim
% levels: for each level in levels crossings will be calculated at scale delta*2^level
% delta: if not 0 then is used as base scale, o/w use median(abs(diff(x)))
% deleteFirst: if 1 then delete first crossing at each level
% isTimeSeries: if 1 then t(n) = n and fname.tim not used
% isIncremental: if 1 then fname.dat contains x(n) - x(n-1), assumes isTimeSeries == 1
% y_eps: jiggle for finding crossings, default value 1e-12
%
% hitting times and points are written to file fname.hit
% subcrossings are written to file fname.sub
%
% Y.Shen 12.2002, O.D.Jones 7.2003, D.A.Rolls 8.2007

% default for y_eps
if nargin < 7, y_eps = 1e-12; end

% hitting times and points
eval(['load ' fname '.dat;']);
eval(['x = ' fname '(:,1);']);
if isTimeSeries
    if isIncremental
        x = [0; cumsum(x)];
    end
    lx = length(x);
    t = 0:lx-1;
else
    lx = length(x);
    eval(['load ' fname '.tim;']);
    eval(['t = ' fname '(:,1);']);
end

if delta == 0
    delta = median(abs(diff(x)));
end

for n = levels
    if n < 10
        output_file = sprintf('%s_0%d.hit', fname, n);
    else
        output_file = sprintf('%s_%d.hit', fname, n);
    end
    fid = fopen(output_file, 'w');
   
    scale = 2^n;
    y = (x - x(1))/delta/scale;
   
    if deleteFirst
        last_hit = 0;
        else
        last_hit = 1;
        end
   
    for i = 2:lx
        if y(i-1) ~= y(i)
            if y(i-1) < y(i)
                x_init = ceil(y(i-1) - y_eps);
                x_final = floor(y(i) + y_eps);
                step = 1;
            else
                x_init = floor(y(i-1) + y_eps);
                x_final = ceil(y(i) - y_eps);
                step = -1;
            end
            for j = x_init:step:x_final
                if j ~= last_hit
                    hit_time = t(i-1) + (j - y(i-1))/(y(i) - y(i-1))*(t(i) - t(i-1));
                    hit_point = j*delta*scale + x(1);
                    fprintf(fid, '%.16f\t%.16f', hit_time, hit_point);
                    fprintf(fid, '\n');
                    last_hit = j;
                end
            end
        end
    end
   
    fclose(fid);
end
   

% subcrossing numbers
levels = levels(2:end);
level = levels(1) - 1;
if level < 10
    fnameXX = sprintf('%s_0%d', fname, level);
else
    fnameXX = sprintf('%s_%d', fname, level);
end
eval(['load ' fnameXX '.hit;']);
eval(['hit0 = ' fnameXX ';']);
   
for level = levels(1:end)
        if level < 10
        fnameXX = sprintf('%s_0%d', fname, level);
        else
        fnameXX = sprintf('%s_%d', fname, level);
        end
        eval(['load ' fnameXX '.hit;']);
    eval(['hit1 = ' fnameXX ';']);
    fid = fopen([fnameXX '.sub'], 'w');
   
    if ~isempty(hit1)
        j0 = 1;
        while hit0(j0,1) ~= hit1(1,1), j0 = j0 + 1; end
        for i = 2:size(hit1,1)
            j1 = j0 + 1;
            while hit0(j1,1) ~= hit1(i,1), j1 = j1 + 1; end
            fprintf(fid, '%d\n', j1 - j0);
            j0 = j1;
        end
        end
   
    fclose(fid);
    hit0 = hit1;
end