关于hurst指数的计算
请教高手以下两种程序求解hurst指数的区别？（本人刚接触matlab，有很多不懂，望赐教）function = estimate_hurst_exponent(data0) % data set
data=data0; % make a local copy
=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 =hurst(d, k)
% Unbiased estimator of the Hurst exponent.
%
% Usage:
% =hurst(d [,k])
%
% INPUTS:
% . d: data
% . k: scales which will be used in the determination.
% . (k may also be of the form or simly 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>
% <a href="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);
if nargin<2
k=;
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=;
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
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
sigma=exp((b1-a12*H)/a11);
if nargout==0
%p=;
=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','vert','bottom','color',)
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',,'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 第一种估算速度比较快
第二种那个文件包应该还有个ar1的文件
MATLAB菜鸟路过~ :handshake
页:
[1]