[Matlab] R/S Hurst exponent function
function [Hal,He,Ht,pval95] = hurst(x,d,fontsize);
%HURST Calculate the Hurst exponent using R/S analysis.
% H = HURST(X) calculates the Hurst exponent of time series X using
% the R/S analysis of Hurst [2], corrected for small sample bias [1,3,4].
% If a vector of increasing natural numbers is given as the second input
% parameter, i.e. HURST(X,D), then it defines the box sizes that the
% sample is divided into (the values in D have to be divisors of the
% length of series X). If D is a scalar (default value D = 50) it is
% treated as the smallest box size that the sample can be divided into.
% In this case the optimal sample size OptN and the vector of divisors
% for this size are automatically computed.
% OptN is defined as the length that possesses the most divisors among
% series shorter than X by no more than 1%. The input series X is
% truncated at the OptN-th value.
% [H,HE,HT] = HURST(X) returns the uncorrected empirical and theoretical
% Hurst exponents.
% [H,HE,HT,PV95] = HURST(X) returns the empirical 95% confidence
% intervals PV95 (see [4]).
%
% If there are no output parameters, the R/S statistics is automatically
% plotted against the divisors on a loglog paper and the results of the
% analysis are displayed in the command window. HURST(X,D,FONTSIZE)
% allows to specify a fontsize different than 14 in the plotted figure.
%
% References:
% [1] A.A.Anis, E.H.Lloyd (1976) The expected value of the adjusted
% rescaled Hurst range of independent normal summands, Biometrica 63,
% 283-298.
% [2] H.E.Hurst (1951) Long-term storage capacity of reservoirs,
% Transactions of the American Society of Civil Engineers 116, 770-808.
% [3] E.E.Peters (1994) Fractal Market Analysis, Wiley.
% [4] R.Weron (2002) Estimating long range dependence: finite sample
% properties and confidence intervals, Physica A 312, 285-299.
% Written by Rafal Weron (2011.09.30).
% Based on functions hurstal.m, hurstcal.m, finddiv.m, findndiv.m
% originally written by Witold Wiland & Rafal Weron (1997.06.30,
% 2001.02.01, 2002.07.27).
if nargin<3,
fontsize = 14;
end
if nargin<2,
d = 50;
end
if max(size(d)) == 1,
% For scalar d set dmin=d and find the 'optimal' vector d
dmin = d;
% Find such a natural number OptN that possesses the largest number of
% divisors among all natural numbers in the interval [0.99*N,N]
N = length(x);
N0 = floor(0.99*N);
dv = zeros(N-N0+1,1);
for i = N0:N,
dv(i-N0+1) = length(divisors(i,dmin));
end
OptN = N0 + find(max(dv)==dv) - 1;
% Use the first OptN values of x for further analysis
x = x(1:OptN);
% Find the divisors of x
d = divisors(OptN,dmin);
else
OptN = length(x);
end
N = length(d);
RSe = zeros(N,1);
ERS = zeros(N,1);
% Calculate empirical R/S
for i=1:N;
RSe(i) = RScalc(x,d(i));
end
% Compute Anis-Lloyd [1] and Peters [3] corrected theoretical E(R/S)
% (see [4] for details)
for i=1:N;
n = d(i);
K = [1:n-1];
ratio = (n-0.5)/n * sum(sqrt((ones(1,n-1)*n-K)./K));
if (n>340)
ERS(i) = ratio/sqrt(0.5*pi*n);
else
ERS(i) = (gamma(0.5*(n-1))*ratio) / (gamma(0.5*n)*sqrt(pi));
end
end
% Calculate the Anis-Lloyd/Peters corrected Hurst exponent
% Compute the Hurst exponent as the slope on a loglog scale
ERSal = sqrt(0.5*pi.*d);
Pal = polyfit(log10(d),log10( RSe - ERS + ERSal ),1);
Hal = Pal(1);
% Calculate the empirical and theoretical Hurst exponents
Pe = polyfit(log10(d),log10(RSe),1);
He = Pe(1);
P = polyfit(log10(d),log10(ERS),1);
Ht = P(1);
% Compute empirical confidence intervals (see [4])
L = log2(OptN);
% R/S-AL (min(divisor)>50) two-sided empirical confidence intervals
pval95 = [0.5-exp(-7.33*log(log(L))+4.21) exp(-7.20*log(log(L))+4.04)+0.5];
C = [ 0.5-exp(-7.35*log(log(L))+4.06) exp(-7.07*log(log(L))+3.75)+0.5 .90];
C = [C; pval95 .95];
C = [C; 0.5-exp(-7.19*log(log(L))+4.34) exp(-7.51*log(log(L))+4.58)+0.5 .99];
% Display and plot results if no output arguments are specified
if nargout < 1,
% Display results
disp('---------------------------------------------------------------')
disp(['R/S-AL using ' num2str(length(d)) ' divisors (' num2str(d(1)) ',...,' num2str(d(length(d))) ...
') for a sample of ' num2str(OptN) ' values'])
disp(['Corrected theoretical Hurst exponent ' num2str(0.5,4)]);
disp(['Corrected empirical Hurst exponent ' num2str(Hal,4)]);
disp(['Theoretical Hurst exponent ' num2str(Ht,4)]);
disp(['Empirical Hurst exponent ' num2str(He,4)]);
disp('---------------------------------------------------------------')
% Display empirical confidence intervals
disp('R/S-AL (min(divisor)>50) two-sided empirical confidence intervals')
disp('--- conf_lo conf_hi level ---------------------------------')
disp(C)
disp('---------------------------------------------------------------')
% Plot R/S
h2 = plot(log10(d),log10(ERSal/(ERS(1)/RSe(1))),'b-');
if fontsize > 10,
set(h2,'linewidth',2);
end;
hold on
h1 = plot(log10(d),log10(RSe-ERS+ERSal),'ro-');
if fontsize > 10,
set(h1,'linewidth',2);
end;
hold off
set(gca,'Box','on','fontsize',fontsize);
xlabel('log_{10}n','fontsize',fontsize);
ylabel('log_{10}R/S','fontsize',fontsize);
legend('Theoretical (R/S)','Empirical (R/S)')
end
function d = divisors(n,n0)
% Find all divisors of the natural number N greater or equal to N0
i = n0:floor(n/2);
d = find((n./i)==floor(n./i))' + n0 - 1;
function rs = RScalc(Z,n)
% Calculate (R/S)_n for given n
m = length(Z)/n;
Y = reshape(Z,n,m);
E = mean(Y);
S = std(Y);
for i=1:m;
Y(:,i) = Y(:,i) - E(i);
end;
Y = cumsum(Y);
% Find the ranges of cummulative series
MM = max(Y) - min(Y);
% Rescale the ranges by the standard deviations
CS = MM./S;
rs = mean(CS);
评论