爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 67132|回复: 60

[程序设计] 有用过重标度极差法计算hurst指数预测趋势的吗?就是R/S法

[复制链接]

新浪微博达人勋

发表于 2015-1-18 21:05:58 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
RSana.m (4.03 KB, 下载次数: 61)
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2015-1-18 22:32:57 | 显示全部楼层

最近看到重标度极差法计算hurst指数可以做趋势预测,所以想研究一下,这是网上找到...

function [logRS,logERS,V]=RSana(x,n,method,q)
% 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
   

这是网上的matlab程序,我有一组50年的数据,怎么求hurst指数来判断未来变化趋势?

RSana.m

4.03 KB, 下载次数: 6, 下载积分: 金钱 -5

matlab

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-3-22 14:35:17 | 显示全部楼层
有一个excel程序,在论坛中搜一下,可以用
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-4-4 11:15:57 | 显示全部楼层
lz,能将你的r/s  excel版本发给我么
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2015-4-4 22:07:01 | 显示全部楼层
启动宏,用那个hurst

RS.xls

56 KB, 下载次数: 142, 下载积分: 金钱 -5

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-4-20 16:15:25 | 显示全部楼层
谢谢楼主分享,没有金钱啊,好东西无法下载(自动回复:请不要使用迅雷等下载工具,点我查看下载帮助
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-5-16 18:07:54 | 显示全部楼层
真是坑啊  下载的不能用
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2015-5-16 22:25:30 | 显示全部楼层
申亚飞 发表于 2015-5-16 18:07
真是坑啊  下载的不能用

5楼的那个,点击启动宏啊
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-5-17 10:13:20 | 显示全部楼层
介绍一下启用宏以后这个软件怎么用呢   有一组时间序列怎么输入?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2015-5-17 11:06:20 | 显示全部楼层
申亚飞 发表于 2015-5-17 10:13
介绍一下启用宏以后这个软件怎么用呢   有一组时间序列怎么输入?

第一列输入你的时间序列,然后点击视图中的宏,再点击执行就可以了
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

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