function [H,sigma]=hurst(d, k)
% Unbiased estimator of the Hurst exponent.
%%

Usage:
% [H,sigma]=hurst(d [,k])
%% INPUTS:
% . d: data
% . k: scales which will be used in the determination.
% . (k may also be of the form [mink maxk] or simly [maxk] which will
% . run faster than explicitly specifying the scales)
%

% INPUTS:
% . H: hurst exponent estimate.
% . sigma: standard dev estimate.
%
% Will make a plot if called with no output arguments.
%% Author: Aslak Grinsted 2007
%% using iterative method described in Koutsoyiannis 2003:
% <a href="http://dx.doi.org/10.1623/hysj.48.1.3.43481"> http://dx.doi.org/10.1623/hysj.48.1.3.43481</a>
% <a href="http://www.itia.ntua.gr/getfile/537/2/2003HSJHurstSuppl.pdf">suppl</a>
% <ahref="http://www.itia.ntua.gr/getfile/537/3/2003HSJHurstPP.pdf">preprint</a>
%% I also really recommend reading this <a href="http://tamino.wordpress.com/2008/06/10/hurst/">blog</a> entry on
Hurst exponents.
% return
if nargin==0
fprintf('No input arguments for hurst. -Loading stockreturns as an
example!...\n')
d=load('stockreturns.mat');d=d.stocks(:,1);
% d=loadproxy('nautadata.txt');d=d(:,2);%hurst(d(:,2))

end
d=d(:);
N=length(d);
i
f nargin<2
k=[1 floor(N/10)];
end
k=k(:);
if length(k)>2
sk=zeros(size(k));
for jj=1:length(k) %slow method but works always....
% ds2=moving(d,k(jj));ds2(isnan(ds2))=[];
% sk(jj,2)=std(ds2(~isnan(ds2)))*k(jj);
ds=filter(ones(k(jj),1),1,d); %moving
ds(1:k(jj)-1)=[];
sk(jj)=std(ds); %dont need to multiply because filter is summing
end
else
if length(k)<2
k=[1 k];
end
if k(1)==1
ds=zeros(size(d,1)+1,1);
else
ds=filter(ones(k(1)-1,1),1,d); %moving
ds(1:k(1)-2)=[];
end
k=(k(1):k(2))';
sk=zeros(size(k));
for jj=1:length(k)
ds(end)=[];
ds=ds+d(k(jj):end); %moving
%sk(jj)=std(ds(1:end-2)+ds(3:end)-ds(2:end-1)*2);
sk(jj)=std(ds);
end
end
lnsk=log(sk);
p=2;
kp=k.^p; %weight from paper
lnk=log(k);
a11=sum(1./kp);
a12=sum(lnk./kp);
%H=polyfit(lnk,lnsk,1); H=H(1); %traditional simplistic estimate
91
APPENDICE C
H=0.5; lastH=inf;
itercount=0;
while abs(H-lastH)>0.00001
lastH=H;
H=min(H,1);
ck=sqrt((N./k-(N./k).^(2*H-1))./(N./k-.5));
lnck=log(ck);
dk=lnk+log(N./k)./(1-(N./k).^(2-2*H));
b1=sum(lnsk./kp)-sum(lnck./kp);
b2=sum(dk.*lnsk./kp)-sum(dk.*lnck./kp);
a21=sum(dk./kp);
a22=sum(dk.*lnk./kp);
H=(a11*b2-a21*b1)/(a11*a22-a21*a12);
itercount=itercount+1;
if itercount>50
error('Hurst.m failed to converge.')
end
end
s
igma=exp((b1-a12*H)/a11);
i
f nargout==0
%p=[H log(sigma)];
[g,a]=ar1(d); vv=a^2/(1-g^2);
g0k=vv * (k.*(1-g^2)-2*g*(1-g.^k))/((1-g)^2) ; %theoretical ar1 ..
eq10 (hurst made easy paper)
fit=ck.*(k.^H)*sigma; %eq.14
loglog(k,sk,'.-',k,fit,'k-',k,sqrt(g0k),'r:')
text(k(round(end*.7)),sk(round(end*.7)),'Empirical','horiz','right','ve
rt','bottom','color',[0 0 1])
text(k(end),fit(end),sprintf('SSS, H=%.2f \\sigma=
%s',H,num2str(sigma,4)),'horiz','right','vert','bottom')
text(k(end),sqrt(g0k(end)),sprintf('AR1, \\gamma=
%.2f',g),'horiz','right','color',[1 0 0],'vert','top')
% text(0.5,0.6,sprintf('H=%.2f \\sigma=
%s',H,num2str(sigma,4)),'horiz','right','units','normalized')
xlabel('k')
ylabel('s_k')
clear H avgsigma sigma
axis tight
end
function H = RS(sequence,isplot)
%%
'RS' estimate the hurst parameter of a given sequence with R/S method.
%%
Inputs:
% sequence: the input sequence for estimate
% isplot: whether display the plot. without a plot if isplot equal to 0
% Outputs:
% H: the estimated hurst coeffeient of the input sequence
% Author: Chu Chen
% Version 1.0, 03/10/2008
% chen-chu@163.com
% i
f nargin == 1
isplot = 0;
end
N = length(sequence);
dlarge = floor(N/5);
dsmall = max(10,log10(N)^2);
D = floor(logspace(log10(dsmall),log10(dlarge),50));
D = unique(D);
n = length(D);
x = zeros(1,n);
y = zeros(1,n);
R = cell(1,n);
S = cell(1,n);
for i = 1:n
d = D(i);
m = floor(N/d);
R{i} = zeros(1,m);
S{i} = zeros(1,m);
matrix_sequence = reshape(sequence(1:d*m),d,m);
Z1 = cumsum(matrix_sequence);
Z2 = cumsum(repmat(mean(matrix_sequence),d,1));
R{i} = (max(Z1-Z2)-min(Z1-Z2));
S{i} = std(matrix_sequence);
if min(R{i})==0 || min(S{i}) ==0
continue;
93
APPENDICE C
end
x(i) = log10(d);
y(i) = mean(log10(R{i}./S{i}));
end
% fit a line with middle part of sequence
index = x~=0;
x = x(index);
y = y(index);
n2 = length(x);
cut_min = ceil(3*n2/10);
cut_max = floor(9*n2/10);
X = x(cut_min:cut_max);
Y = y(cut_min:cut_max);
p1 = polyfit(X,Y,1);
Yfit = polyval(p1,X);
H = (Yfit(end)-Yfit(1))/(X(end)-X(1));
f isplot ~= 0
figure,hold on;
bound = ceil(log10(N));
axis([0 bound 0 0.75*bound]);
temp = (1:n).*index;
index = temp(index);
for i = 1:n2
plot(x(i),log10(R{index(i)}./S{index(i)}),'b.');
end
x = linspace(0,bound,10);
y1 = 0.5*x;
y2 = x;
h1 = plot(x,y1,'b--','LineWidth',2);
h2 = plot(x,y2,'b-.','LineWidth',2);
plot(X,Yfit,'r-','LineWidth',3);
legend([h1,h2],'slope 1/2','slope 1',4)
xlabel('log10(blocks of size m)'),ylabel('log10(R/S)'),title('R/S
Method');
end

Successivo: Modelli Arch e Garch

Sommario: Indice