马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
x
以前用程序计算过一个序列的Hurst指数，现在matlab里重新算，出现下列错误，请问各位高手如何让解决？请高手帮我运行下！以前运行可以的，电脑坏了，重装了matlab,可是应该可以啊？纠结啊！！！
function [logRS,logERS,V]=RSana(x,n,method,q)
%Syntax: [logRS,logERS,V]=RSana(x,n,method,q)
%____________________________________________
%
% Performs R/S analysis on a time series.
%
% logRS is the log(R/S).
% logERS is the Expectation of log(R/S).
% V is the V statistic.
% x is the time series.
% n is the vector with the sub-periods.
% method can take one of the following values
% 'Hurst' for the Hurst-Mandelbrot variation.
% 'Lo' for the Lo variation.
% 'MW' for the Moody-Wu variation.
% 'Parzen' for the Parzen variation.
% q can be either
% a (non-negative) integer.
% 'auto' for the Lo's suggested value.
%
%
% References:
%
% Peters E (1991): Chaos and Order in the Capital Markets. Willey
%
% Peters E (1996): Fractal Market Analysis. Wiley
%
% Lo A (1991): Long term memory in stock market prices. Econometrica
% 59: 1279-1313
%
% Moody J, Wu L (1996): Improved estimates for Rescaled Range and Hurst
% exponents. Neural Networks in Financial Engineering, eds. Refenes A-P
% Abu-Mustafa Y, Moody J, Weigend A: 537-553, Word Scientific
%
% Hauser M (1997): Semiparametric and nonparametric testing for long
% memory: A Monte Carlo study. Empirical Economics 22: 247-271
%
%
% Alexandros Leontitsis
% Department of Education
% University of Ioannina
% 45110 - Dourouti
% Ioannina
% Greece
%
% University e-mail: me00743@cc.uoi.gr
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% 1 Jan 2004.
if nargin<1 | isempty(x)==1
error('You should provide a time series.');
else
% x must be a vector
if min(size(x))>1
error('Invalid time series.');
end
x=x(:);
% N is the time series length
N=length(x);
end
if nargin<2 | isempty(n)==1
n=1;
else
% n must be either a scalar or a vector
if min(size(n))>1
error('n must be either a scalar or a vector.');
end
% n must be integer
if n-round(n)~=0
error('n must be integer.');
end
% n must be positive
if n<=0
error('n must be positive.');
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 must be a scalar
if sum(size(q))>2
error('q must be scalar.');
end
% q must be integer
if q-round(q)~=0
error('q must be integer.');
end
% q must be positive
if q<0
error('q must be positive.');
end
end
end
for i=1:length(n)
% Calculate the sub-periods
a=floor(N/n(i));
% Make the sub-periods matrix
X=reshape(x(1:a*n(i)),n(i),a);
% Estimate the mean of each sub-period
ave=mean(X);
% Remove the mean from each sub-period
cumdev=X-ones(n(i),1)*ave;
% Estimate the cumulative deviation from the mean
cumdev=cumsum(cumdev);
% Estimate the standard deviation
switch method
case 'Hurst'
% Hurst-Mandelbrot variation
stdev=std(X);
case 'Lo'
% Lo variation
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 variation
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 variation
if mod(q,2)~=0
error('For the "Parzen" variation q must be dived by 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('You should provide another value for "method".');
end
% Estiamte the rescaled range
rs=(max(cumdev)-min(cumdev))./stdev;
clear stdev
% Take the logarithm of the mean R/S
logRS(i,1)=log10(mean(rs));
if nargout>1
% Initial calculations fro the log(E(R/S))
j=1:n(i)-1;
s=sqrt((n(i)-j)./j);
s=sum(s);
% The estimation of log(E(R/S))
logERS(i,1)=log10(s/sqrt(n(i)*pi/2));
% Other estimations of 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
% Estimate V
V(i,1)=mean(rs)/sqrt(n(i));
end
end
调用
q=0;x=10.386
10.703
11.055
11.413
11.753
12.059
12.33
12.58
12.832
13.111
13.438
13.818
14.245
14.693
15.133
15.531
15.866
16.13
16.333
16.5
16.665
16.864
17.119
17.44
17.817
18.224
18.627
18.993
19.295
19.524
19.689
19.811
19.923
20.056
20.231
20.458
20.73
21.026
21.322
21.593
21.824
22.012
22.17
22.321
22.49
22.701
22.965
23.283
23.642
24.018
24.388
24.733
25.044
25.324
25.589
25.861
26.162
26.507
26.901
27.337
27.796
28.256
28.694
29.096
29.461
29.796
30.121
30.456
30.821
31.226
31.671
32.142
32.621
33.084
33.516
33.906
34.259
34.588
34.913
35.255
35.629
36.042
36.488
36.953
37.417
37.86
38.271
38.644
38.986
39.313
39.641
39.988
40.364
40.771
41.199
41.635
42.063
42.469
42.847
43.2
43.54
43.884
44.247
44.642
45.073
45.533
46.01
46.488
46.951
47.388
47.797
48.187
48.571
48.964
49.383
49.835
50.321
50.834
51.36
51.884
52.395
52.887
53.362
53.829
54.3
54.787
55.298
55.835
56.39
56.952
57.505
58.037
58.542
59.018
59.475
59.924
60.382
60.861
61.366
61.898
62.449
63.006
63.558
64.093
64.609
65.109
65.6
66.095
66.604
67.133
67.682
68.244
68.808
69.362
69.896
70.405
70.891
71.363
71.833
72.317
72.823
73.357
73.917
74.496
75.081
75.663
76.233
76.79
77.337
77.884
78.441
79.018
79.617
80.238
80.873
81.511
82.139
82.75
83.338
83.907
84.465
85.024
85.595
86.186
86.801
87.435
88.081
88.727
89.364
89.986
90.592
91.189
91.784
92.388
93.009
93.651
;n=1:196;method='Hurst';
[logRS,logERS,V]=RSana(x,n,method,q)
错误如下
Warning: Divide by zero.
> In mean at 29
In RSana at 201
Warning: Divide by zero.
> In mean at 29
In RSana at 181
Warning: Divide by zero.
> In mean at 29
In RSana at 201
Warning: Divide by zero.
> In mean at 29
In RSana at 181
Warning: Divide by zero.
> In mean at 29
In RSana at 201
[ 本帖最后由 parklyn 于 2010-4-26 17:07 编辑 ] |
发表于 2010-4-26 16:53
提升卡
变色卡
显身卡
发表于 2010-4-26 16:57
楼主
