eda_staty.m

计算时间序列的Hurst系数有许多方法

   %j2 = length(nj);
   %j1 =  newchoosej1(regu,nj, muj, 0,j2 );
   %j1= max(subj1);
   j1=j1full;
   j2=j2full;
   [alphaest,cfCest,cfest,Cest,Q,Valpha,VcfC,CoValphacfC,Vcf,CoValphacf,unsafe,yj,varj,aest] = regrescomp(regu,nj, muj, j1,j2,0);
   fullalpha = alphaest;
   fullstdalpha = sqrt(Valpha);
   fullH = (alphaest + 1)/2;
   fullstdH = sqrt(Valpha)/2;  
   fullcf = cfest;
   fullstdcf = sqrt(Vcf);
   CIalpha = fullalpha;
   CIalpha = max(CIalpha,0.01);
   CIalpha = min(CIalpha,0.9);
   cgam = cfest./( 2*((2*pi).^CIalpha) .* gamma(CIalpha) .* sin(pi*(1-CIalpha)/2) )   ;
   LRDasympstdfull = sqrt(cgam)./( lensub.^((1-CIalpha)/2) .* sqrt((1+CIalpha).*CIalpha/2) )  ;
   
   if wantH   % set output variables correctly
      al = H;
      stdal = stdH;
      fullal = fullH;
      fullstdal = fullstdH;
   else
      al = alpha;
      stdal = stdalpha;
      fullal = fullalpha;
      fullstdal = fullstdalpha;
   end 

   if (printout)
     figure(20)
   end
end
   


if (lookmean)  
   m = mean(series');
   stdm = std(series');
% output for Balaton,
seuil*stdm(1)/sqrt(lensub)
seuil*LRDasympstd(1)
   if (printout)
     nplotnum = nplotnum+1;
     subplot(nplotnum)
     plot(1:nsub,m,'k-')
     hold
     plot(1:nsub,mean(m).*line,'k-')
     for n = 1:nsub
%seuil*stdm(n)/sqrt(lensub)
       plot([n n],[m(n)-seuil*stdm(n)/sqrt(lensub) m(n)+seuil*stdm(n)/sqrt(lensub) ],'b-') ;   % add SRD confidence intervals
       if ( H(n)> 1/2 & H(n) <1)
         plot([n+.05 n+.05],[m(n)-seuil*LRDasympstd(n) m(n)+seuil*LRDasympstd(n) ],'r-');% add asymptotic LRD CI's
       else
         plot([n+.05 n+.05],[m(n)-seuil*stdm(n)/sqrt(lensub) m(n)+seuil*stdm(n)/sqrt(lensub) ],'r-');  % SRD CI's if Hest<1/2.
       end
     end
     V = axis;
     axis([ 0.5 nsub+0.5 V(3) V(4) ]); 
     %xlabel('Subseries number')
     %xlabel('Block number')
     %ylabel('Mean rates, megabits/s')
     ylabel('Means')
    % title('Means of the subseries. The horizontal lines give the overall value and confidence intervals.')
     title('Means over the blocks. The horizontal line gives the overall mean.')
     grid

     %  Add CI's for entire trace on the left of mean plot.
     meanfull = mean(m) ;
     stdfull  = std(x)  ;
     plot([.7 .7],[meanfull-seuil*stdfull/sqrt(len) meanfull+seuil*stdfull/sqrt(len) ],'b-');  % SRD to compare
     if ( fullH> 1/2 & fullH <1)
       plot([.75 .75],[meanfull-seuil*LRDasympstdfull meanfull+seuil*LRDasympstdfull ],'r-');   % LRD
     %  plot([.8,1:nsub],meanfull+seuil*LRDasympstdfull.*CIline,'r-')  %  top line of CI 
     %  plot([.8,1:nsub],meanfull-seuil*LRDasympstdfull.*CIline,'r-')  %  bottom line of CI
     else
       plot([.7 .7],[meanfull-seuil*stdfull/sqrt(len) meanfull+seuil*stdfull/sqrt(len) ],'b-');% SRD if H<1/2 or>1
     end
     hold off
   end
end
if (lookvar)
   v = (std(series')).^2;
   vfull = std(x)^2;
   %  if X~N(mu,v), then the sample var (biased) is ~ v/n * ChiSq_n-1, with var = 2v^2 *(n-1)/n^2
   %  assume then that for large n the sample var is Gaussian with var 2v^2/n
   stdasympSRD = sqrt(2/lensub)*v ;   % using v to estimate the real v for each subseries
   if (printout)
     nplotnum = nplotnum+1;
     subplot(nplotnum)
     %plot(log2(1:nsub),log2(v),'-*')
     plot(1:nsub,v,'k-')
     hold
     for n = 1:nsub
       plot([n n],[v(n)-seuil*stdasympSRD(n) v(n)+seuil*stdasympSRD(n)],'b-') ;   %SRD to compare
     end
     plot(1:nsub,mean(v)*line,'--')    %  add average of v's of the subseries as a dashed line
     
     %  Add CI's for entire trace on the left of variance plot
     stdasympSRD = sqrt(2/len)*vfull ;
     plot([.75 .75],[vfull-seuil*stdasympSRD vfull+seuil*stdasympSRD],'b-') ;
     %plot([.8,1:nsub],vfull+seuil*stdasympSRD.*CIline,'b-')  %  top line of CI
     %plot([.8,1:nsub],vfull-seuil*stdasympSRD.*CIline,'b-')  %  bottom line of CI
     plot(1:nsub,vfull*line,'k-')  %  add overall variance as a solid line
     
     V = axis;
     axis([ 0.5 nsub+0.5 V(3) V(4) ]); 
     hold off
     %xlabel('Subseries number')
     ylabel('Variances')
     %title('Variance of the subseries. The solid (dashed) horizontal line gives the overall (average) value.')
     grid
   end
end
if (lookH)
   if (printout)
     nplotnum = nplotnum+1;
     subplot(nplotnum)
  
     if (wantH)
       plot(1:nsub,H,'k-')
       HL = H-seuil*stdH;
       HR = H+seuil*stdH;
       fullHL = fullH-seuil*fullstdH;
       fullHR = fullH+seuil*fullstdH;
       hold on
       for n = 1:nsub
         plot([n n],[HL(n) HR(n)],'r-') ;   % add confidence intervals
       end
       plot([.75 .75],[fullHL fullHL],'r-') 
       plot(1:nsub,mean(H)*line,'--')    %  add average of H's of the subseries as a dashed line
     else
       plot(1:nsub,alpha,'k-')
       alphaL = alpha-seuil*stdalpha;
       alphaR = alpha+seuil*stdalpha;
       fullalphaL = fullalpha-seuil*fullstdalpha;
       fullalphaR = fullalpha+seuil*fullstdalpha; 
       hold on
       for n = 1:nsub
         plot([n n],[alphaL(n) alphaR(n)],'r-') ;   % add confidence intervals
       end
       plot([.75 .75],[fullalphaL fullalphaR],'r-')
       plot(1:nsub,mean(alpha)*line,'--')   
     end
     
     if (wantH)
       fprintf('For the full series, (j1,j2)= (%d,%d),  H = %4.2f,  cf = %7.4f,  Q= %5.4f \n\n',j1,j2,fullH,cf,Q)
       plot(1:nsub,fullH*line,'k-')
       %plot([.8,1:nsub],fullH+seuil*fullstdH.*CIline,'r-')  %  top line of CI
       %plot([.8,1:nsub],fullH-seuil*fullstdH.*CIline,'r-')  %  bottom line of CI
       V = axis;
       axis([ 0.5 nsub+0.5 min(HL)-0.2*abs(min(HL)-mean(H)) max(HR)+0.2*abs(min(HL)-mean(H)) ]);
       %xlabel('Subseries number')
       ylabel('H')
       title('Scaling parameter of the subseries. The solid (dashed) horizontal line gives the overall (average) value.')
     else
       fprintf('For the full series, (j1,j2)= (%d,%d),  alpha = %4.2f,   cf = %7.4f,  Q= %5.4f \n\n',j1,j2,fullalpha,fullcf,Q)
       plot(1:nsub,fullalpha*line,'k-')
       %plot([.8,1:nsub],fullalpha+seuil*fullstdalpha.*CIline,'r-')  %  top line of CI
       %plot([.8,1:nsub],fullalpha-seuil*fullstdalpha.*CIline,'r-')  %  bottom line of CI
       V = axis;
       axis([ 0.5 nsub+0.5 min(alphaL)-0.2*abs(min(alphaL)-mean(alpha)) max(alphaR)+0.2*abs(min(alphaL)-mean(alpha)) ]);
       ylabel('\alpha')
       title('Scaling parameter of the subseries. The solid (dashed) horizontal line gives the overall (average) value.')
     end
     grid
   end

   %%%%%%  Constancy of alpha test
   if nsub>1 
    siglevel = 0.05;
    sigma = stdalpha(1);  % all the std's are the same for same size blocks
    mean(alpha);
    %std(alpha)
    %std(alpha)^2
    %%%%% Calculate for general test
    V = ( (std(alpha))^2*(nsub-1)/nsub )*nsub/sigma^2 ;   % statistic = S^2 * m/sig^2 
    %% the crit level is passed back, it can be compared against siglevel in the calling program
    crit_level = 1 - chi2_CDF( V, nsub-1); % test probability corresponding exactly to the data
 
    %%%%% Calculate for two-level test form each value of m1, in [1,nsub-1]
    crit_level2 = ones(1,nsub-1);
    m1 = max(1,floor(cutpoints/4));   % trick for A4 section 4.3, only calculate the one needed : m/4
    %m1 = max(1,floor(cutpoints/2));   % trick for A4 section 4.2, only calculate the one needed : m/2
    % for m1 = 1:(nsub-1)
      xbar = mean(alpha(1:m1));
      ybar = mean(alpha((m1+1):nsub));
      V2(m1) = abs(xbar-ybar)/sigma;
      scaledsigma(m1) = sqrt( nsub/m1/(nsub-m1) );
      crit_level2(m1) = 2*gauss_CDF( -V2(m1), 0, scaledsigma(m1) );   % gauss_CDF can't take vector q
    %end
    %V2;

    if printout
      if wantH        % use this for legend placement
         textposition = max(HR);
      else
         textposition = max(alphaR); 
      end 
      C = chi2_CDFinv(1-siglevel,nsub-1);           % boundary of critical region 
      if  V < C
        fprintf('\n***** General constancy test for alpha cannot be rejected at significance level %4.3f \n', siglevel)
        Handle=text(nsub*0.8-1+0.1,textposition,['\bf Not Rejected, (',num2str(100*(crit_level),3), '% sig)']) ;
        %Handle=text(nsub*0.8-1+0.1,textposition,'\bf Not Rejected') ;
        set(Handle,'fontsize',12);
      else
        fprintf('\n***** General constancy test for alpha was rejected at significance level %4.3f \n',siglevel)
        Handle=text(nsub*0.8-1+0.1,textposition,['\bf Rejected, (',num2str(100*(crit_level),3), '% sig)'] ) ;
        %       Handle=text(nsub*0.8-1+0.1,textposition,'\bf Rejected') ;
        set(Handle,'fontsize',12); 
      end
      fprintf('      The test statistic was %4.2f corresponding to a critical level of %4.3f and probability of %4.3f \n',V,crit_level,1-crit_level)
      fprintf('      The critical region was to the right of %4.3f  \n', C)
    
      %% two-level output
      C2 = gauss_CDFinv(1-siglevel/2, 0, scaledsigma );    % critical regions get bigger with m1
      fprintf('\n      The two-level tests, beginning at m1 = 1,  gave critical levels of: \n')
      crit_level2
      fprintf('      Corresponding rejections ( 1=reject ) were:\n')
      rejections = (V2 >  C2);
      rejections

      hold off
    end  

   end  % constancy test

   %ans=input('Should I launch LDestimate on the trace? (hit return for yes)  ');
   %if  isempty(ans)
      %LDestimate(x,regu,1,j2,1)
   %end
end


if (lookcf)
   if (printout)
     nplotnum = nplotnum+1;
     subplot(nplotnum)
     plot(1:nsub,cf,'k-')
     sig_level = 5 ;  %set confidence level
     m = log( (cf.^4)./((stdcf.^2)+cf.*cf) )/2;
     v = log( (stdcf.^2)./cf./cf + 1 );
     cfL = exp(gauss_CDFinv(  sig_level/2/100,m,sqrt(v)));
     cfR = exp(gauss_CDFinv(1-sig_level/2/100,m,sqrt(v)));

     m = log( (fullcf^4)/((fullstdcf^2)+fullcf*fullcf) )/2;
     v = log( (fullstdcf^2)/fullcf/fullcf + 1 );
     fullcfL = exp(gauss_CDFinv(  sig_level/2/100,m,sqrt(v)));
     fullcfR = exp(gauss_CDFinv(1-sig_level/2/100,m,sqrt(v)));

     hold on
     for n = 1:nsub
       plot([n n],[cfL(n) cfR(n) ],'r-') ;   % add confidence intervals
     end
     plot([.75 .75],[fullcfL fullcfR],'r-')
     plot(1:nsub,mean(cf)*line,'--')    %  add average of cf's of the subseries as a dashed line

     plot(1:nsub,fullcf*line,'k-')
     %plot([.8,1:nsub],fullcf+seuil*fullstdcf.*CIline,'r-')  %  top line of CI
     %plot([.8,1:nsub],fullcf-seuil*fullstdcf.*CIline,'r-')  %  bottom line of CI
     V = axis;
     axis([ 0.5 nsub+0.5 min(cfL)-1.6*abs(min(cfL)-mean(cf)), max(cfR)+0.8*abs(min(cfR)-mean(cf)) ]);
     %xlabel('Subseries number')
     ylabel('c_f')
     title('Second scaling parameter of the subseries. The solid (dashed) horizontal line gives the overall (average) value.')
     hold off
     grid
     fprintf('For the full series, (j1,j2)= (%d,%d),  alpha = %4.2f,  cf = %7.4f,   Q= %5.4f \n\n',j1,j2,fullalpha, fullcf,Q)
   end
end

if (printout)
  xlabel('Block number')
  fprintf('\n*******************************************************************************************\n  \n')
end