RS分形算法matlab代码实现

happy

RS分形算法matlab代码实现,既利用RS方法计算时间序列的分形hurst指数

1. function [logRS,logERS,V]=RSana(x,n,method,q)
   
2. %Syntax: [logRS,logERS,V]=RSana(x,n,method,q)
   
3. %____________________________________________
   
4. %
   
5. % Performs R/S analysis on a time series.
   
6. %
   
7. % logRS is the log(R/S).
   
8. % logERS is the Expectation of log(R/S).
   
9. % V is the V statistic.
   
10. % x is the time series.
    
11. % n is the vector with the sub-periods.
    
12. % method can take one of the following values
    
13. %  'Hurst' for the Hurst-Mandelbrot variation.
    
14. %  'Lo' for the Lo variation.
    
15. %  'MW' for the Moody-Wu variation.
    
16. %  'Parzen' for the Parzen variation.
    
17. % q can be either
    
18. %  a (non-negative) integer.
    
19. %  'auto' for the Lo's suggested value.
    
20. %
    
21. %
    
22. % References:
    
23. %
    
24. % Peters E (1991): Chaos and Order in the Capital Markets. Willey
    
25. %
    
26. % Peters E (1996): Fractal Market Analysis. Wiley
    
27. %
    
28. % Lo A (1991): Long term memory in stock market prices. Econometrica
    
29. % 59: 1279-1313
    
30. %
    
31. % Moody J, Wu L (1996): Improved estimates for Rescaled Range and Hurst
    
32. % exponents. Neural Networks in Financial Engineering, eds. Refenes A-P
    
33. % Abu-Mustafa Y, Moody J, Weigend A: 537-553, Word Scientific
    
34. %
    
35. % Hauser M (1997): Semiparametric and nonparametric testing for long
    
36. % memory: A Monte Carlo study. Empirical Economics 22: 247-271
    
37. %
    
38. %
    
39. % Alexandros Leontitsis
    
40. % Department of Education
    
41. % University of Ioannina
    
42. % 45110 - Dourouti
    
43. % Ioannina
    
44. % Greece
    
45. %
    
46. % University e-mail: me00743@cc.uoi.gr
    
47. % Lifetime e-mail: leoaleq@yahoo.com
    
48. % Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
    
49. %
    
50. % 1 Jan 2004.
    
51. 
    
52. if nargin<1 | isempty(x)==1
    
53.    error('You should provide a time series.');
    
54. else
    
55.    % x must be a vector
    
56.    if min(size(x))>1
    
57.       error('Invalid time series.');
    
58.    end
    
59.    x=x(:);
    
60.    % N is the time series length
    
61.    N=length(x);
    
62. end
    
63. 
    
64. if nargin<2 | isempty(n)==1
    
65.    n=1;
    
66. else
    
67.    % n must be either a scalar or a vector
    
68.    if min(size(n))>1
    
69.       error('n must be either a scalar or a vector.');
    
70.    end
    
71.    % n must be integer
    
72.    if n-round(n)~=0
    
73.        error('n must be integer.');
    
74.    end
    
75.    % n must be positive
    
76.    if n<=0
    
77.       error('n must be positive.');
    
78.    end
    
79. end
    
80. 
    
81. if nargin<4 | isempty(q)==1
    
82.    q=0;
    
83. else
    
84.     if q=='auto'
    
85.         t=autocorr(x,1);
    
86.         t=t(2);
    
87.         q=((3*N/2)^(1/3))*(2*t/(1-t^2))^(2/3);
    
88.     else
    
89.         % q must be a scalar
    
90.         if sum(size(q))>2
    
91.             error('q must be scalar.');
    
92.         end
    
93.         % q must be integer
    
94.         if q-round(q)~=0
    
95.             error('q must be integer.');
    
96.         end
    
97.         % q must be positive
    
98.         if q<0
    
99.             error('q must be positive.');
    
100.         end
     
101.     end
     
102. end
     
103. 
     
104. 
     
105. for i=1:length(n)
     
106.    
     
107.     % Calculate the sub-periods
     
108.     a=floor(N/n(i));
     
109.    
     
110.     % Make the sub-periods matrix
     
111.     X=reshape(x(1:a*n(i)),n(i),a);
     
112.    
     
113.     % Estimate the mean of each sub-period
     
114.     ave=mean(X);
     
115.    
     
116.     % Remove the mean from each sub-period
     
117.     cumdev=X-ones(n(i),1)*ave;
     
118.    
     
119.     % Estimate the cumulative deviation from the mean
     
120.     cumdev=cumsum(cumdev);
     
121.    
     
122.     % Estimate the standard deviation
     
123.     switch method
     
124.     case 'Hurst'
     
125.         % Hurst-Mandelbrot variation
     
126.         stdev=std(X);
     
127.     case 'Lo'
     
128.         % Lo variation
     
129.         for j=1:a
     
130.             sq=0;
     
131.             for k=0:q
     
132.                 v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
     
133.                 if k>0
     
134.                     sq=sq+(1-k/(q+1))*v(k+1);
     
135.                 end
     
136.             end
     
137.             stdev(j)=sqrt(v(1)+2*sq);
     
138.         end
     
139.     case 'MW'
     
140.         % Moody-Wu variation
     
141.         for j=1:a
     
142.             sq1=0;
     
143.             sq2=0;
     
144.             for k=0:q
     
145.                 v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
     
146.                 if k>0
     
147.                     sq1=sq1+(1-k/(q+1))*(n(i)-k)/n(i)/n(i);
     
148.                     sq2=sq2+(1-k/(q+1))*v(k+1);
     
149.                 end
     
150.             end
     
151.             stdev(j)=sqrt((1+2*sq1)*v(1)+2*sq2);
     
152.         end
     
153.     case 'Parzen'
     
154.         % Parzen variation
     
155.         if mod(q,2)~=0
     
156.             error('For the "Parzen" variation q must be dived by 2.');
     
157.         end
     
158.         for j=1:a
     
159.             sq1=0;
     
160.             sq2=0;
     
161.             for k=0:q
     
162.                 v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
     
163.                 if k>0 & k<=q/2
     
164.                     sq1=sq1+(1-6*(k/q)^2+6*(k/q)^3)*v(k+1);
     
165.                 elseif k>0 & k>q/2
     
166.                     sq2=sq2+(1-(k/q)^3)*v(k+1);
     
167.                 end
     
168.             end
     
169.             stdev(j)=sqrt(v(1)+2*sq1+2*sq2);
     
170.         end
     
171.     otherwise
     
172.         error('You should provide another value for "method".');
     
173.     end
     
174.    
     
175.     % Estiamte the rescaled range
     
176.     rs=(max(cumdev)-min(cumdev))./stdev;
     
177.    
     
178.     clear stdev
     
179.    
     
180.     % Take the logarithm of the mean R/S
     
181.     logRS(i,1)=log10(mean(rs));
     
182.    
     
183.     if nargout>1
     
184.         
     
185.         % Initial calculations fro the log(E(R/S))
     
186.         j=1:n(i)-1;
     
187.         s=sqrt((n(i)-j)./j);
     
188.         s=sum(s);
     
189.         
     
190.         % The estimation of log(E(R/S))
     
191.         logERS(i,1)=log10(s/sqrt(n(i)*pi/2));
     
192.         
     
193.         % Other estimations of log(E(R/S))
     
194.         %logERS(i,1)=log10((n(i)-0.5)/n(i)*s/sqrt(n(i)*pi/2));
     
195.         %logERS(i,1)=log10(sqrt(n(i)*pi/2));
     
196.         
     
197.     end
     
198.    
     
199.     if nargout>2
     
200.         % Estimate V
     
201.         V(i,1)=mean(rs)/sqrt(n(i));
     
202.     end
     
203. 
     
204. end

复制代码

linapeng

你好，请问这个代码输出hurst外，其他几个是甚么呢

suffer

原帖由 linapeng 于 2006-11-27 12:52 发表
你好，请问这个代码输出hurst外，其他几个是甚么呢

% logRS is the log(R/S).
% logERS is the Expectation of log(R/S).
% V is the V statistic.

自己看注释吧

pheigenbau

这么强的程序啊?

zjamy

分形hurst指数
是什么意思呢？
就是分形维数吗？

suffer

原帖由 pheigenbau 于 2006-11-30 16:53 发表
这么强的程序啊?

你指什么？

suffer

原帖由 zjamy 于 2006-12-4 21:49 发表
分形hurst指数
是什么意思呢？
就是分形维数吗？

H.E.HURST（赫斯特）是英国水文学家。以他命名的HURST指数，被广泛用于资本市场的混沌分形分析。除了埃德加.E.彼得斯的两本专著外，近几年也发表了一些论文。
一个具有赫斯特统计特性的系统，不需要通常概率统计学的独立随机事件假设。它反映的是一长串相互联系事件的结果。今天发生的事将影响未来，过去的事也会影响现在。这正是我们分析资本市场所需要的理论和方法。传统的概率统计学，对此是难办到的。
HURST指数（H）有三个不同类型：
1、H=0.5，标志着所研究的序列是一个随机序列，即过去的增量与未来的增量不相关。这是通常概率统计学的研究对象；
2、0.5<H<1.0，标志着所研究的序列是一个持久性序列，即过去的增量与未来的增量正相关。序列有长程相关性；
3、0<H<0.5，标志着所研究的序列是一个反持久性序列，即过去的增量与未来的增量负相关，序列有突变跳跃逆转性。
根据赫斯特的研究，自然界的很多自然现象，H大于0.5。埃德加.E.彼得斯的两本专著，对国外资本市场进行了系统分析，证实了许多市场指数的H也大于0.5；近几年国内发表了一些论文，同样验证了沪深市场指数的H也大于0.5。这种市场特征，被称为是有偏随机游动市场，也即市场具有混沌分形特征。
对于反持久性序列，埃德加.E.彼得斯指出，它在经济金融学中虽然很重要，但人们发现的反持久性序列却很少。因此，在两本专著中论述不多。
HURST指数的经典计算方法，是R/S分析法，即重标极差分析法。用此法计算HURST指数，不仅计算量大，且方法繁杂。目前所见论文，一般都是针对少数代表性指数，且多半是用月（周）数据分析的。对于个股的HURST指数计算，尚未见到。在现有的几个股软（飞狐或分析家）上直接实现，虽有可能，也较困难。因此，除了经典计算方法外，寻求一种简单些的计算方法也是有必要的。

为了便于对沪深市场的所有股票进行分析，本人采用了一种简化方法，计算出沪深市场所有股票的HURST指数。对这种方法只进行了一般性的误差分析，未与R/S分析法进行全面比较。因此，计算结果和初步认识仅供参考。
一、指数类的HURST指数
1、上证：0.568，沪A：0.571，沪B：0.583；
2、深成指：0.533，深综指0.515
初步认识：就大盘而论，沪深市场的日线周期数据序列是持久性序列，且沪市比深市有偏性更强，沪B比沪A有偏性更强。
二、A股市场所有股票的HURST指数（1260支）
1、H大于等于0.55的有36支，H大于等于0.52的有96支，H大于0.50的有137支
2、H小于等于0.45的有36支，H小于等于0.48的有1016支
初步认识：A股市场的绝大多数股票，其日线周期数据序列是反持久性序列，序列有极强的突变逆转性。相反，持久性序列仅有1/10左右。这是现有研究中无人提及的。
三、B股市场所有股票的HURST指数（109支）
H大于0.50的有69支，H小于0.50的有31支
初步认识：B股市场的股票比A股市场有偏性强，但仍有3/10的B股，其日线周期数据序列是反持久性序列。
说明：新股的日线周期数据序列太短，按理是不能计算HURST指数的。统计时未剔除新股。

以上引自MACD股市技术分析俱乐部

wsmrzw

这是上文中提到的简化算法吗

lavender25

这个程序如何输入数据，运行这个.m文件后出现
??? Error using ==> RSana123
You should provide a time series.
还有就是输出的logRS is the log(R/S)应该不是最后的hurst指数吧

furoo

为什么输入数据，出来这个结果？程序有误？
??? Error: File: C:\MATLAB7\work\Rsana.m Line: 1 Column: 21
Unbalanced or misused parentheses or brackets.

furoo

这个程序只能输出一个log(rs)，并且好像就是hurst指数吧，因为就一个值，其他hurs期望指数和V统计值输不出来，加入代码输出总是显示错误

furoo

原帖由 furoo 于 2009-2-17 20:02 发表

                               
登录/注册后可看大图

为什么输入数据，出来这个结果？程序有误？
??? Error: File: C:\MATLAB7\work\Rsana.m Line: 1 Column: 21
Unbalanced or misused parentheses or brackets.

这个问题已经解决了，是我第一行注释没搞好

evita1986

好东西，学习！

parklyn

不是的 分形维数和hurst指数是完全不同的两个概念 而且此程序出来的结果不是hurst指数 将数据拟合后得到直线的斜率是hurst指数

cr810317

楼主人太太太太好了，感谢分享:loveliness: