|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
用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);
[logRS,logERS,V]=RSana(x,n,'Hurst',1);
运行后得到的logRS和logERS带有NaN和-Inf,这是什么意思?
函数代码如下:
- function [logRS,logERS,V]=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
复制代码
|
|