分形参数计算程序分享计算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的默认值.ifnargin<1|isempty(x)==1error('你应该给出一个时间序列.');else%x必须是变量ifmin(size(x))>1error('时间序列无效.');endx=x(:);%N是时间序列的长度N=length(x);endifnargin<2|isempty(n)==1n=1;else%n必须是一个变化的标量或矢量ifmin(size(n))>1error('n必须是一个变化的标量或矢量.');end%n必须是个整数ifn-round(n)~=0error('n必须是个整数.');end%n必须是确定ifn<=0error('n必须是确定.');endendifnargin<4|isempty(q)==1q=0;elseifq=='auto't=autocorr(x,1);t=t(2);q=((3*N/2)^(1/3))*(2*t/(1-t^2))^(2/3);else%q必须是标量ifsum(size(q))>2error('q必须是标量.');end%q必须是整数ifq-round(q)~=0error('q必须是整数.');end%q必须是确定ifq<0error('q必须是确定.');endendendfori=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);%估算这个标准偏差switchmethodcase'Hurst'%Hurst-Mandelbrot参数stdev=std(X);case'Lo'%Lo参数forj=1:asq=0;fork=0:qv(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);ifk>0sq=sq+(1-k/(q+1))*v(k+1);endendstdev(j)=sqrt(v(1)+2*sq);endcase'MW'%Moody-Wu参数forj=1:asq1=0;sq2=0;fork=0:qv(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);ifk>0sq1=sq1+(1-k/(q+1))*(n(i)-k)/n(i)/n(i);sq2=sq2+(1-k/(q+1))*v(k+1);endendstdev(j)=sqrt((1+2*sq1)*v(1)+2*sq2);endcase'Parzen'%Parzen参数ifmod(q,2)~=0error('在"Parzen"参数中q必须是2.');endforj=1:asq1=0;sq2=0;fork=0:qv(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1);ifk>0&k<=q/2sq1=sq1+(1-6*(k/q)^2+6*(k/q)^3)*v(k+1);elseifk>0&k>q/2sq2=sq2+(1-(k/q)^3)*v(k+1);endendstdev(j)=sqrt(v(1)+2*sq1+2*sq2);endotherwiseerror('你应该付给"method"另一个值.');end%估算(R(t)/s(t))rs=(max(cumdev)-min(cumdev))./stdev;clearstdev%取这个平均值R/S的对数logRS(i,1)=log10(mean(rs));ifnargout>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));endifnargout>2%估算VV(i,1)=mean(rs)/sqrt(n(i));endend计算lyapunovCODE:functionlambda_1=lyapunov_wolf(data,N,m,tau,P)%该函数用来计算时间序列的最大Lyapunov指数--Wolf方法%m:嵌入维数%tau:时间延迟%data:时间序列%N:时间序列长度%P:时间序列的平均周期,选择演化相点距当前点的位置差，即若当前相点为I，则演化相点只能在|I－J|>P的相点中搜寻%lambda_1:返回最大lyapunov指数值min_point=1;%&&要求最少搜索到的点数MAX_CISHU=5;%&&最大增加搜索范围次数%FLYINGHAWK%求最大、最小和平均相点距离max_d=0;%最大相点距离min_d=1.0e+100;%最小相点距离avg_dd=0;Y=reconstitution(data,N,m,tau);%相空间重构M=N-(m-1)*tau;%重构相空间中相点的个数fori=1:(M-1)forj=i+1:Md=0;fork=1:md=d+(Y(k,i)-Y(k,j))*(Y(k,i)-Y(k,j));endd=sqrt(d);ifmax_ddmin_d=d;endavg_dd=avg_dd+d;endendavg_d=2*avg_dd/(M*(M-1));%平均相点距离dlt_eps=(avg_d-min_d)*0.02;%若在min_eps～max_eps中找不到演化相点时，对max_eps的放宽幅度min_eps=min_d+dlt_eps/2;%演化相点与当前相点距离的最小限max_eps=min_d+2*dlt_eps;%&&演化相点与当前相点距离的最大限%从P+1～M-1个相点中找与第一个...
