【急】RS分形算法计算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);
[logRS,logERS,V]=RSana(x,n,'Hurst',1);


运行后得到的logRS和logERS带有NaN和-Inf，这是什么意思？

           




函数代码如下：

1. function [logRS,logERS,V]=RSana(x,n,method,q)
   
2. % 用 R/S 方法分析间序列
   
3. % logRS 是 log(R/S).
   
4. % logERS 是 log(R/S)期望.
   
5. % V 是统计量.
   
6. % x 是时间序列.
   
7. % n 是这个数列的子集.
   
8. % method 可以取下列值
   
9.   %  'Hurst' 为了Hurst-Mandelbrot变量
   
10.   %  'Lo' 是Lo变量.
    
11.   %  'MW' 是Moody-Wu变量.
    
12.   %  'Parzen' 是Parzen变量.
    
13. % q 可以是任意值
    
14. %  a 是非0整数.
    
15. %  'auto' 是 Lo的默认值.
    
16. 
    
17. if nargin<1 | isempty(x)==1
    
18.    error('你应该给出一个时间序列.');
    
19. else
    
20.    % x 必须是变量
    
21.    if min(size(x))>1
    
22.       error('时间序列无效.');
    
23.    end
    
24.    x=x(:);
    
25.    % N 是时间序列的长度
    
26.    N=length(x);
    
27. end
    
28. 
    
29. if nargin<2 | isempty(n)==1
    
30.    n=1;
    
31. else
    
32.    % n 必须是一个变化的标量或矢量
    
33.    if min(size(n))>1
    
34.       error('n 必须是一个变化的标量或矢量.');
    
35.    end
    
36.    % n 必须是个整数
    
37.    if n-round(n)~=0
    
38.        error('n 必须是个整数.');
    
39.    end
    
40.    % n 必须是确定
    
41.    if n<=0
    
42.       error('n 必须是确定.');
    
43.    end
    
44. end
    
45. 
    
46. if nargin<4 | isempty(q)==1
    
47.    q=0;
    
48. else
    
49.     if q=='auto'
    
50.         t=autocorr(x,1);
    
51.         t=t(2);
    
52.         q=((3*N/2)^(1/3))*(2*t/(1-t^2))^(2/3);
    
53.     else
    
54.         % q 必须是标量
    
55.         if sum(size(q))>2
    
56.             error('q 必须是标量.');
    
57.         end
    
58.         % q 必须是整数
    
59.         if q-round(q)~=0
    
60.             error('q 必须是整数.');
    
61.         end
    
62.         % q 必须是确定
    
63.         if q<0
    
64.             error('q 必须是确定.');
    
65.         end
    
66.     end
    
67. end
    
68. 
    
69. 
    
70. for i=1:length(n)
    
71.    
    
72.     % 计算这个子序列
    
73.     a=floor(N/n(i));
    
74.    
    
75.     % 仓健这个子序列的矩阵
    
76.     X=reshape(x(1:a*n(i)),n(i),a);
    
77.    
    
78.     % 估算这个子序列的平均值
    
79.     ave=mean(X);
    
80.    
    
81.     % 给这个序列的每一个值除以平均值
    
82.     cumdev=X-ones(n(i),1)*ave;
    
83.    
    
84.     % 估算累计离差
    
85.     cumdev=cumsum(cumdev);
    
86.    
    
87.     % 估算这个标准偏差
    
88.     switch method
    
89.     case 'Hurst'
    
90.         % Hurst-Mandelbrot 参数
    
91.         stdev=std(X);
    
92.     case 'Lo'
    
93.         % Lo 参数
    
94.         for j=1:a
    
95.             sq=0;
    
96.             for k=0:q
    
97.                 v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
    
98.                 if k>0
    
99.                     sq=sq+(1-k/(q+1))*v(k+1);
    
100.                 end
     
101.             end
     
102.             stdev(j)=sqrt(v(1)+2*sq);
     
103.         end
     
104.     case 'MW'
     
105.         % Moody-Wu 参数
     
106.         for j=1:a
     
107.             sq1=0;
     
108.             sq2=0;
     
109.             for k=0:q
     
110.                 v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
     
111.                 if k>0
     
112.                     sq1=sq1+(1-k/(q+1))*(n(i)-k)/n(i)/n(i);
     
113.                     sq2=sq2+(1-k/(q+1))*v(k+1);
     
114.                 end
     
115.             end
     
116.             stdev(j)=sqrt((1+2*sq1)*v(1)+2*sq2);
     
117.         end
     
118.     case 'Parzen'
     
119.         % Parzen 参数
     
120.         if mod(q,2)~=0
     
121.             error('在"Parzen" 参数中q 必须是2.');
     
122.         end
     
123.         for j=1:a
     
124.             sq1=0;
     
125.             sq2=0;
     
126.             for k=0:q
     
127.                 v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
     
128.                 if k>0 & k<=q/2
     
129.                     sq1=sq1+(1-6*(k/q)^2+6*(k/q)^3)*v(k+1);
     
130.                 elseif k>0 & k>q/2
     
131.                     sq2=sq2+(1-(k/q)^3)*v(k+1);
     
132.                 end
     
133.             end
     
134.             stdev(j)=sqrt(v(1)+2*sq1+2*sq2);
     
135.         end
     
136.     otherwise
     
137.         error('你应该付给 "method"另一个值.');
     
138.     end
     
139.    
     
140.     % 估算(R(t)/s(t))
     
141.     rs=(max(cumdev)-min(cumdev))./stdev;
     
142.    
     
143.     clear stdev
     
144.    
     
145.     % 取这个平均值 R/S的对数
     
146. 
     
147.     logRS(i,1)=log10(mean(rs));
     
148.    
     
149.     if nargout>1
     
150.         
     
151.         % 开始计算log(E(R/S))
     
152.         j=1:n(i)-1;
     
153.         s=sqrt((n(i)-j)./j);
     
154.         s=sum(s);
     
155.         
     
156.         % 估算log(E(R/S))
     
157.         logERS(i,1)=log10(s/sqrt(n(i)*pi/2));
     
158.         
     
159.         %其它估算log(E(R/S))
     
160.         %logERS(i,1)=log10((n(i)-0.5)/n(i)*s/sqrt(n(i)*pi/2));
     
161.         %logERS(i,1)=log10(sqrt(n(i)*pi/2));
     
162.         
     
163.     end
     
164.    
     
165.     if nargout>2
     
166.         % 估算 V
     
167.         V(i,1)=mean(rs)/sqrt(n(i));
     
168.     end
     
169. 
     
170. end
     

复制代码