分形参数计算程序分享.doc

分形参数计算程序分享计算hurst指数 CODE: 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 nargin1 | isempty(x)==1? ?error(你应该给出一个时间序列.);else? ?% x 必须是变量? ?if min(size(x))1? ?? ?error(时间序列无效.);? ?end? ?x=x(:);? ?% N 是时间序列的长度? ?N=length(x);endif nargin2 | 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 必须是确定.);? ?endendif nargin4 | 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 q0? ?? ?? ?? ?error(q 必须是确定.);? ?? ???end? ? endendfor 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 k0? ?? ?? ?? ?? ?? ???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 k0? ?? ??