欢迎来到天天文库
浏览记录
ID:34442301
大小:48.67 KB
页数:10页
时间:2019-03-06
《分形参数计算程序分享》由会员上传分享,免费在线阅读,更多相关内容在工程资料-天天文库。
1、分形参数计算程序分享计算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
2、isempty(x)==1 e
3、rror('你应该给出一个时间序列.');else %x必须是变量 ifmin(size(x))>1 error('时间序列无效.'); end x=x(:); %N是时间序列的长度 N=length(x);endifnargin<2
4、isempty(n)==1 n=1;else %n必须是一个变化的标量或矢量 ifmin(size(n))>1 error('n必须是一个变化的标量或矢量.'); end %n必须是个整数 ifn-round(n)~=0 error('n必须是个整数.'); end %n必须是确定 ifn<=0 error
5、('n必须是确定.'); endendifnargin<4
6、isempty(q)==1 q=0;else ifq=='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))>2 error('q必须是标量.'); end %q必须是整数 ifq-round(q)~=0 error('q必须是整数.'); end %q必须
7、是确定 ifq<0 error('q必须是确定.'); end endendfori=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); %估算这个标准偏差 switchmethod case'Hur
8、st' %Hurst-Mandelbrot参数 stdev=std(X); case'Lo' %Lo参数 forj=1:a sq=0; fork=0:q v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1); ifk>0 sq=sq+(1-k/(q+1))*v(k+1); end end stdev(j)=sqrt(v(1)+2*sq); end case'M
9、W' %Moody-Wu参数 forj=1:a sq1=0; sq2=0; fork=0:q v(k+1)=sum(X(k+1:n(i),j)'*X(1:n(i)-k,j))/(n(i)-1); ifk>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*sq
10、1)*v(1)+2*sq2); end case'Parzen' %Parzen参数 ifmod(q,2)~=0 error('在"Parzen"参数中q必须是2.'); end forj=1:a sq1=0; sq2=0; fork=0:q v(k+1)=sum(X(k+1:n(i)
此文档下载收益归作者所有