function [logRS,logERS,V]=RSana(x,n,method,q)
% 用 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
   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(x,1);
        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)-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 參數
        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 參數
        if mod(q,2)~=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)-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('你應該付給 "method"另一個值.');
    end
    
    % 估算(R(t)/s(t))
    rs=(max(cumdev)-min(cumdev))./stdev;
    
    clear stdev
    
    % 取這個平均值 R/S的對數

    logRS(i,1)=log10(mean(rs));
    
    if nargout>1
        
        % 開始計算log(E(R/S))
        j=1:n(i)-1;
        s=sqrt((n(i)-j)./j);
        s=sum(s);
        
        % 估算log(E(R/S))
        logERS(i,1)=log10(s/sqrt(n(i)*pi/2));
        
        %其它估算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
        % 估算 V
        V(i,1)=mean(rs)/sqrt(n(i));
    end

end