关于R/S分析程序用法

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

看到有人用下面的源代码，但为何有时候求出的值是负数，请各位帮忙，我是新接触最好能稍微具体一点讲解，谢谢。

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
       x(i)=log10(i);
        hold on
        plot(x(i),V(i,1),'-')
        plot(x(i),logRS(i,1),'g')
        plot(x(i),logERS(i,1),'r')
end

请问程序完整吗，为什么我得不到logERS和V统计量呢？:L

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享1 收藏2 回帖

关键词：分析程序 Mandelbrot otherwise function nargout 程序

%Syntax: [logRS,logERS,V]=RSana(x,n,method,q)
% x is the time series.
% n is the vector with the sub-periods.
% method can take one of the following values
% 'Hurst' for the Hurst-Mandelbrot variation.
%  'Lo' for the Lo variation.
%  'MW' for the Moody-Wu variation.
%  'Parzen' for the Parzen variation.
% q can be either
% a (non-negative) integer.
%  'auto' for the Lo's suggested value.
x=randn(3000,1);
x=cumsum(x);
[logRS,logERS,V]=RSana(x,100,'Hurst',1)

logRS =

    1.5696


logERS =

    1.0557


V =

    3.7117

>> x=randn(3000,1);
>> x=cumsum(x);
>> [logRS,logERS,V]=RSana(x,100,'Hurst',1)

logRS =

    5.2480


logERS =

    3.5068


V =

    3.8001

怎么差异怎么大？
还有其它参数 比如 ‘Lo’, 'MW','Parzen' 没有问题吗， 我用'Hurst'一般没出现问题，但其它选其它参数就出现负数，不知道问题在哪？

谢谢。

yonghu123 发表于 2012-4-1 23:34
>> x=randn(3000,1);
>> x=cumsum(x);
>> [logRS,logERS,V]=RSana(x,100,'Hurst',1)

数据由乱数产生,我的结果自然跟你的不同
四种method:'Hurst','Lo','MW','Parzen'都OK的

谢谢了。

多谢~~受教了·

你好，我最近在做这个，但是不知道那个程序里面的n如何设定，能否告知？加QQ指导更好，564637620~~

epoh 发表于 2012-4-2 08:36
数据由乱数产生,我的结果自然跟你的不同
四种method:'Hurst','Lo','MW','Parzen'都OK的

你好，能不能讲一下就是这个里面的n的设置是如何设定啊？是乱设置的么？我看楼主的X是3000行，设值n=100,这个是可以任意定的么？

epoh 发表于 2012-4-2 08:36
数据由乱数产生,我的结果自然跟你的不同
四种method:'Hurst','Lo','MW','Parzen'都OK的

请问MW和Parzen具体是怎么求的？有相关文献介绍一下计算的过程吗？

yonghu123 发表于 2012-4-1 23:34
>> x=randn(3000,1);
>> x=cumsum(x);
>> [logRS,logERS,V]=RSana(x,100,'Hurst',1)

请问一下，x=cumsum(x)，为什么要进行累加呢？