声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1405|回复: 0

[综合] 【急】用R/S分形算法计算hurst指数遇到问题?

[复制链接]
发表于 2013-5-19 15:12 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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,这是什么意思?


捕获.PNG     捕获2.PNG


函数代码如下:
  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. if nargin<1 | isempty(x)==1
  17.    error('你应该给出一个时间序列.');
  18. else
  19.    % x 必须是变量
  20.    if min(size(x))>1
  21.       error('时间序列无效.');
  22.    end
  23.    x=x(:);
  24.    % N 是时间序列的长度
  25.    N=length(x);
  26. end

  27. if nargin<2 | isempty(n)==1
  28.    n=1;
  29. else
  30.    % n 必须是一个变化的标量或矢量
  31.    if min(size(n))>1
  32.       error('n 必须是一个变化的标量或矢量.');
  33.    end
  34.    % n 必须是个整数
  35.    if n-round(n)~=0
  36.        error('n 必须是个整数.');
  37.    end
  38.    % n 必须是确定
  39.    if n<=0
  40.       error('n 必须是确定.');
  41.    end
  42. end

  43. if nargin<4 | isempty(q)==1
  44.    q=0;
  45. else
  46.     if q=='auto'
  47.         t=autocorr(x,1);
  48.         t=t(2);
  49.         q=((3*N/2)^(1/3))*(2*t/(1-t^2))^(2/3);
  50.     else
  51.         % q 必须是标量
  52.         if sum(size(q))>2
  53.             error('q 必须是标量.');
  54.         end
  55.         % q 必须是整数
  56.         if q-round(q)~=0
  57.             error('q 必须是整数.');
  58.         end
  59.         % q 必须是确定
  60.         if q<0
  61.             error('q 必须是确定.');
  62.         end
  63.     end
  64. end


  65. for i=1:length(n)
  66.    
  67.     % 计算这个子序列
  68.     a=floor(N/n(i));
  69.    
  70.     % 仓健这个子序列的矩阵
  71.     X=reshape(x(1:a*n(i)),n(i),a);
  72.    
  73.     % 估算这个子序列的平均值
  74.     ave=mean(X);
  75.    
  76.     % 给这个序列的每一个值除以平均值
  77.     cumdev=X-ones(n(i),1)*ave;
  78.    
  79.     % 估算累计离差
  80.     cumdev=cumsum(cumdev);
  81.    
  82.     % 估算这个标准偏差
  83.     switch method
  84.     case 'Hurst'
  85.         % Hurst-Mandelbrot 参数
  86.         stdev=std(X);
  87.     case 'Lo'
  88.         % Lo 参数
  89.         for j=1:a
  90.             sq=0;
  91.             for k=0:q
  92.                 v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
  93.                 if k>0
  94.                     sq=sq+(1-k/(q+1))*v(k+1);
  95.                 end
  96.             end
  97.             stdev(j)=sqrt(v(1)+2*sq);
  98.         end
  99.     case 'MW'
  100.         % Moody-Wu 参数
  101.         for j=1:a
  102.             sq1=0;
  103.             sq2=0;
  104.             for k=0:q
  105.                 v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
  106.                 if k>0
  107.                     sq1=sq1+(1-k/(q+1))*(n(i)-k)/n(i)/n(i);
  108.                     sq2=sq2+(1-k/(q+1))*v(k+1);
  109.                 end
  110.             end
  111.             stdev(j)=sqrt((1+2*sq1)*v(1)+2*sq2);
  112.         end
  113.     case 'Parzen'
  114.         % Parzen 参数
  115.         if mod(q,2)~=0
  116.             error('在"Parzen" 参数中q 必须是2.');
  117.         end
  118.         for j=1:a
  119.             sq1=0;
  120.             sq2=0;
  121.             for k=0:q
  122.                 v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);
  123.                 if k>0 & k<=q/2
  124.                     sq1=sq1+(1-6*(k/q)^2+6*(k/q)^3)*v(k+1);
  125.                 elseif k>0 & k>q/2
  126.                     sq2=sq2+(1-(k/q)^3)*v(k+1);
  127.                 end
  128.             end
  129.             stdev(j)=sqrt(v(1)+2*sq1+2*sq2);
  130.         end
  131.     otherwise
  132.         error('你应该付给 "method"另一个值.');
  133.     end
  134.    
  135.     % 估算(R(t)/s(t))
  136.     rs=(max(cumdev)-min(cumdev))./stdev;
  137.    
  138.     clear stdev
  139.    
  140.     % 取这个平均值 R/S的对数

  141.     logRS(i,1)=log10(mean(rs));
  142.    
  143.     if nargout>1
  144.         
  145.         % 开始计算log(E(R/S))
  146.         j=1:n(i)-1;
  147.         s=sqrt((n(i)-j)./j);
  148.         s=sum(s);
  149.         
  150.         % 估算log(E(R/S))
  151.         logERS(i,1)=log10(s/sqrt(n(i)*pi/2));
  152.         
  153.         %其它估算log(E(R/S))
  154.         %logERS(i,1)=log10((n(i)-0.5)/n(i)*s/sqrt(n(i)*pi/2));
  155.         %logERS(i,1)=log10(sqrt(n(i)*pi/2));
  156.         
  157.     end
  158.    
  159.     if nargout>2
  160.         % 估算 V
  161.         V(i,1)=mean(rs)/sqrt(n(i));
  162.     end

  163. end
复制代码


本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-11 01:00 , Processed in 0.066435 second(s), 23 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表