【急】用R/S分形算法计算hurst指数遇到问题？
用RS分形算法计算Hurst指数：（代码来自http://forum.vibunion.com/thread-32137-1-1.html）运行：x=xlsread('C:\Users\yearn\Desktop\1\S170Q2.00.xlsx','A1:A2000');
n=1:length(x);
=RSana(x,n,'Hurst',1);
运行后得到的logRS和logERS带有NaN和-Inf，这是什么意思？
函数代码如下：
function =RSana(x,n,method,q)
% 用 R/S 方法分析间序列
% logRS 是 log(R/S).
% logERS 是 log(R/S)期望.
% V 是统计量.
% x 是时间序列.
% n 是这个数列的子集.
% method 可以取下列值
%'Hurst' 为了Hurst-Mandelbrot变量
%'Lo' 是Lo变量.
%'MW' 是Moody-Wu变量.
%'Parzen' 是Parzen变量.
% q 可以是任意值
%a 是非0整数.
%'auto' 是 Lo的默认值.
if nargin<1 | isempty(x)==1
error('你应该给出一个时间序列.');
else
% x 必须是变量
if min(size(x))>1
error('时间序列无效.');
end
x=x(:);
% N 是时间序列的长度
N=length(x);
end
if nargin<2 | isempty(n)==1
n=1;
else
% n 必须是一个变化的标量或矢量
if min(size(n))>1
error('n 必须是一个变化的标量或矢量.');
end
% n 必须是个整数
if n-round(n)~=0
error('n 必须是个整数.');
end
% n 必须是确定
if n<=0
error('n 必须是确定.');
end
end
if nargin<4 | isempty(q)==1
q=0;
else
if q=='auto'
t=autocorr(x,1);
t=t(2);
q=((3*N/2)^(1/3))*(2*t/(1-t^2))^(2/3);
else
% q 必须是标量
if sum(size(q))>2
error('q 必须是标量.');
end
% q 必须是整数
if q-round(q)~=0
error('q 必须是整数.');
end
% q 必须是确定
if q<0
error('q 必须是确定.');
end
end
end
for i=1:length(n)
% 计算这个子序列
a=floor(N/n(i));
% 仓健这个子序列的矩阵
X=reshape(x(1:a*n(i)),n(i),a);
% 估算这个子序列的平均值
ave=mean(X);
% 给这个序列的每一个值除以平均值
cumdev=X-ones(n(i),1)*ave;
% 估算累计离差
cumdev=cumsum(cumdev);
% 估算这个标准偏差
switch method
case 'Hurst'
% Hurst-Mandelbrot 参数
stdev=std(X);
case 'Lo'
% Lo 参数
for j=1:a
sq=0;
for k=0:q
v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
if k>0
sq=sq+(1-k/(q+1))*v(k+1);
end
end
stdev(j)=sqrt(v(1)+2*sq);
end
case 'MW'
% Moody-Wu 参数
for j=1:a
sq1=0;
sq2=0;
for k=0:q
v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
if k>0
sq1=sq1+(1-k/(q+1))*(n(i)-k)/n(i)/n(i);
sq2=sq2+(1-k/(q+1))*v(k+1);
end
end
stdev(j)=sqrt((1+2*sq1)*v(1)+2*sq2);
end
case 'Parzen'
% Parzen 参数
if mod(q,2)~=0
error('在"Parzen" 参数中q 必须是2.');
end
for j=1:a
sq1=0;
sq2=0;
for k=0:q
v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
if k>0 & k<=q/2
sq1=sq1+(1-6*(k/q)^2+6*(k/q)^3)*v(k+1);
elseif k>0 & k>q/2
sq2=sq2+(1-(k/q)^3)*v(k+1);
end
end
stdev(j)=sqrt(v(1)+2*sq1+2*sq2);
end
otherwise
error('你应该付给 "method"另一个值.');
end
% 估算(R(t)/s(t))
rs=(max(cumdev)-min(cumdev))./stdev;
clear stdev
% 取这个平均值 R/S的对数
logRS(i,1)=log10(mean(rs));
if nargout>1
% 开始计算log(E(R/S))
j=1:n(i)-1;
s=sqrt((n(i)-j)./j);
s=sum(s);
% 估算log(E(R/S))
logERS(i,1)=log10(s/sqrt(n(i)*pi/2));
%其它估算log(E(R/S))
%logERS(i,1)=log10((n(i)-0.5)/n(i)*s/sqrt(n(i)*pi/2));
%logERS(i,1)=log10(sqrt(n(i)*pi/2));
end
if nargout>2
% 估算 V
V(i,1)=mean(rs)/sqrt(n(i));
end
end
页:
[1]