RS分析--重标极差分析法程序

function [logRSlogERSV]=RSana（xnmethodq）
% 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（x1）;
        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;