重标极差分析法R/S分析

R/S分析法也称重标极差分析法，同时引入了一个统计量Hurst指数。Hurst指数常用于分析时间序列的分形特征和长期记忆过程，目前在时间序列变化趋势的持续性或反持续性强度判断方面得到广泛引用。

function [logRSlogERSVH]=RSana（xnmethodq）
% 用 R/S 方法分析间序列
% logRS 是 log（R/S）.
% logERS 是 log（R/S）期望.
% V 是统计量.
% x 是时间序列.
% n 是这个数列的子集.
% method 可以取下列值
%  ‘Hurst‘ 为了Hurst-Mandelbrot变量
%  ‘Lo‘ 是Lo变量.
%  ‘MW‘ 是Moody-Wu变量.
%  ‘Parzen‘ 是Parzen变量.
% q 可以是任意值
%  a 是非0整数.
%  ‘auto‘ 是 Lo的默认值.

if nargin<1 | isempty（x）==1     %nargin为函数输入变量数
   error（‘你应该给出一个时间序列.‘）;
else
   
   % x 必须是变量
   if min（size（x））>1
      error（‘时间序列无效.‘）;
   end
   x=x（:）;

   % N 是时间序列的长度
   N=length（x）;
end



if nargin<2 | isempty（n）==1
   n=1;
else
   % n 必须是一个变化的标量或矢量
   if min（size（n））>1
      error（‘n 必须是一个变化的标量或矢量.‘）;
   end
   
   % n 必须是个整数
   if n-round（n）~=0
       error（‘n 必须是个整数.‘）;
   end

   % n 必须是确定
   if n<=0
      error（‘n 必须是确定.‘）;
   end
end



if nargin<4 | isempty（q）==1
   q=0;
else
    if q==‘auto‘
        t=autocorr（x1）;
        t=t（2）;
        q=（（3*N/2）^（1/3））*（2*t/（1-t^2））^（2/3）;
    else
        % q 必须是标量
        if sum（size（q））>2
            error（‘q 必须是标量.‘）;
        end
        % q 必须是整数
        if q-round（q）~=0
            error（‘q 必须是整数.‘）;
        end
        % q 必须是确定
        if q<0
            error（‘q 必须是确定.‘）;
        end
    end
end



for i=1:length（n）
    
    % 计算这个子序列
    a=floor（N/n（i））;
    
    % 创建这个子序列的矩阵
    X=reshape（x（1:a*n（i））n（i）a）;
    
    % 估算这个子序列的平均值
    ave=mean（X）;
    
    % 给这个序列的每一个值除以平均值 
    cumdev=X-ones（n（i）1）*ave;
    
    % 估算累计离差
    cumdev=cumsum（cumdev）;
    
    % 估算这个标准偏差
    switch method
    case ‘Hurst‘
        % Hurst-Mandelbrot 参数
        stdev=std（X）;
    case ‘Lo‘
        % Lo 参数
        for j=1:a
            sq=0;
            for k=0:q
                v（k+1）=sum（X（k+1:n（i）j）‘*X（1:n（i）-kj））/（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 参数
        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）-kj））/（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 参数
        if mod（q2）~=0
            error（‘在“Parzen“ 参数中q 必须是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）-kj））/（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（‘你应该付给 “method“另一个值.‘）;
    end
    
   

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件         758  2009-05-22 23:02  rs_analysis.m
     文件        4075  2009-05-22 23:12  RSana.m
     文件        4572  2009-05-22 23:12  run.asv
     文件        4547  2009-05-22 23:17  runmovie.m
     文件        4644  2009-05-22 23:46  runnews.m