eda_staty.m

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                       %
%   eda_staty.m                                         %
%                                                       %              
%        D. Veitch   P.Abry                             %
%                                                       %
%   25/1/98                                             %
%   DV, Lyon  14 oct                                    %
%   DV, Melbourne, 2 May 2000                           %
%   DV, Melbourne, 24 Dec 2000                          %
%   DV, Melb, 10 Sep 2004:  robustified j1,j2 choice    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%    This function supports an experimental analysis of stationarity. 
%    The series is cut into subseries and tests performed on each for:
%        - the mean   , with confidence intervals: LRD and IID
%        - the variance
%        - H   and the stationarity test  (alpha or H are possible, see the internal variable "wantH" 
%                                          which selects which one)
%        - cf
%    
%    Other choices can easily be added!! 
%
%             eg:  eda_staty(delay45,1,1,1,1,1,20,  2   ,2,6,2,6, 1,1);   
%                  Examine delay45  cut into twenty equal size 
%                  subseries, considering mean, var, and H with N=2, and plotting the series
%                  Test is for the regime [j1,j2]=[2,6]
%                  eda_staty(fgn8     ,1,1,0,1,1, 5,  2   ,2,7,2,7  , 1, 1);  % with initialisation j1*=2 !
%                  eda_staty(data     ,1,1,0,0,0,  6,  2   ,6,13,6,13, 0, 1);  % demo for CI's of mean
%                  eda_staty(pOct_work',0,1,0,1,0, 4,  2   ,6,14,6,14, 0, 1);  % acceptance of constancy of alpha     
%
%
%--- Routines called Directly :
%
%   % wtspec, newchoosej1,  regrescomp, initDWT_discrete, chi2_CDF, chi2_CDFinv, gauss_CDF, gauss_CDFinv
%
%  Input:   x:    the series to be analysed (is converted to a row vector)
%           looksig :   logical parameter, 1 means the signal will be plotted
%           lookmean:   logical parameter, 1 means the mean will be  analysed, else should be zero.
%           lookvar :   logical parameter, 1 means the variance   "   "
%           lookH   :   logical parameter, 1 means the exponent (H or alpha)   "   "
%           lookcf  :   logical parameter, 1 means the cf S       "   " 
%           cutpoints:  value or vector describing how the series should be partitioned.
%                         - if a value=nsub, then nsub subseries are created of equal size.
%                         - else a vector of specific values is taken (*** not currently supported***) 
%           regu:  (regularity) =  number of vanishing moments (of the Daubechies wavelet).
%           j1slice:   the j1 value for each block
%           j2slice:   the j2 value for each block
%           j1full:   the j1 value for the entire series
%           j2full:   the j2 value for the entire series
%           discrete_init: 1:  perform the special MRA initialisation for intrinsically discrete series
%                         else: assume data already initialised (ie is already the approximation sequence)
%           printout:   1: -text output and the multiplot graph, one row per statistic selected.
%                          -if LD's are calculated, then outputted in a plot, one for each block, Q in title.
%                    else: nothing.
%
%  Output:  m:   the vector of means of each subseries
%           v:   the vector of variances of each subseries
%           al:   the vector of scaling exponents (H or alpha) values of each subseries 
%           cf:    the vector of cf values of each subseries
%           stdal: the corresponding vector of standard deviations of the estimates (theoretical)
%           stdcf: the corresponding vector of standard deviations of the estimates (theoretical)
%           fullal:     the exponent value (H or alpha) over the entire series
%           fullcf:     the cf value over the entire series
%           fullstdal:  the corresponding  standard deviation of the estimate (theoretical)
%           fullstdcf: the corresponding  standard deviation of the estimate (theoretical)
%           crit_level:  the 'complementary' quantile (expressed as a probability) of the constantcy test statistic, ie probability that the data is at least this bad.
%           crit_level2: a vector of values for the two-level test with m1 = 1..nsub/2  .
%
%  A set of subplots are plotted on figure 20.  First the series, if selected, followed by 
%  a subplot for each quantity selected 
%
%  Call: [m,v,al,cf,stdal,stdcf,fullal,fullcf,fullstdcf,fullstdal,crit_level,crit_level2] = eda_staty(x, looksig, lookmean, lookvar, lookH, lookcf, cutpoints , regu, j1slice, j2slice, j1full, j2full, discrete_init, printout)
%

function [m,v,al,cf,stdal,stdcf,fullal,fullcf,fullstdcf,fullstdal,crit_level,crit_level2] = eda_staty(x, looksig, lookmean, lookvar, lookH, lookcf, cutpoints , regu, j1slice, j2slice, j1full, j2full, discrete_init, printout)

%%%%%%%%%% internal variables
wantH = 0;   % print and output according to H or alpha,   alpha = 2H - 1


%%%%%%%%%%% correct input
if  (lookH | lookmean | lookcf)  % if want to look at H, make sure regu is set.
   regu=max(regu,1);
end
dim = size(x);
if dim(1) > 1
  x = x';
end

%%%%%%%%%%% init
x = x';
len = length(x);
m = 0;
v = 0;
al = 0;
stdal = 0;
fullal = 0;
fullstdal = 0;
crit_level = 0;


%%%%%%%%%%%  determine subseries of x
if (length(cutpoints)==1)
  %cutpoints = len*(1:cutpoints)/cutpoints
  %nsub = length(cutpoints)
  nsub  = cutpoints;
  lensub = floor(len/nsub);
  %cutpoints
  %endpts =  (1:nsub)*lensub ;    % output indices of endpoints of each block
  %endpts(endpts==4424)  % check special value for A4
else
  nsub = length(cutpoints) + 1;
  lensub = diff([1 cutpoints len]) % *** after this unequal lengths are not supported! but could be done
end 
   

%%%%%%%%%%  create subseries as rows of matrix 'series'
series = ones(nsub,lensub);
for n = 1:nsub
  series(n,:) = x( (1+(n-1)*lensub):(n*lensub) )';
end

series(:,1:min(12,lensub));

if ( (lookmean | lookvar | lookH | lookcf ) & printout)
fprintf('\n*******************************************************************************************\n\n')
  fprintf('%d subseries of length %d will be examined\n\n',nsub,lensub)
end


%%%%%%%%%%%  For each quantity, calculate the value for each subseries
seuil = 1.9599 ;
num_plot = ( looksig + lookmean + lookvar + lookH + lookcf );
nplotnum = (num_plot)*100 + 10;
if (printout)
  figure(20)
  clf
end
if (looksig & printout)
   nplotnum = nplotnum+1;
   subplot(nplotnum)
   skip = ceil(len/5000);
   plot(1:skip:len,x(1:skip:len),'k-');   % don't plot all the points of huge timeseries
   axis tight
   %title('Timeseries:  fGn of length 2^{18},  \alpha = 0.6,  c_f = 100')
   title('Timeseries')
   %title('Timeseries, pOct aggregated in 10ms bins')
end

line   = ones(1,nsub);     
CIline = ones(1,nsub+1);

if (lookH | lookcf | lookmean)   %  need H values For confidence intervals for mean as well as for H itself
   alpha = [];
   H = [];
   cf = [];
   subj1 = [];
   subQ = [];

   if (printout)
     fprintf('\n\n');
     figure(21)
     clf
     minlower = 100000000;    %  variables to fit vertical axis in LD plots
     minlowerfull = minlower;
     maxupper= -minlower;
     maxupperfull = maxupper;
   end

   for n = 1:nsub
      %  Initialize the MRA based recursive method of calculating the DWT coefficients
      if  discrete_init==1    % use the special initialisation for intrinsically discrete series
         filterlength = 0;    %  choose the automatic length selection algorithm, or set here if desired
         [appro,kfirst,klast] = initDWT_discrete(series(n,:),regu,filterlength,0);      % no output
         if printout
            fprintf('** Using initialization for discrete series, filterlength = %d\n',lensub+kfirst-klast)
         end
         %n = klast-kfirst+1;
         [muj,nj]=wtspec(appro,regu,fix( log2(klast-kfirst+1) ));    % use length of appro, after filtering
         if printout 
           figure(21);
         end 
      else
         appro = series(n,:);
         if printout
            fprintf('** Taking the given data as the initial approximation sequence\n')
         end
         [muj,nj]=wtspec(series(n,:),regu,fix( log2(lensub) ));
      end

   %  j2 = length(nj);
      j2 = min(j2slice,length(nj)) ;
      if (j2<2 )
        fprintf('Subseries too short to evaluate H! aborting.  \n');
        %j1 = j1slice;
        %H = 1/2.*line;                                                                  
        return;
      else
%        j1 =  newchoosej1(regu,nj, muj, 0,j2 );  % no printout, just the LRD case.
        if (j1slice>=j2)
          fprintf('Subseries too short to evaluate H! please change (j1,j2) selection. For now, trying j1-1 \n');
          j1 = max(j2-1,1);
        else
          j1 = j1slice;
        end
        subj1 = [ subj1 j1 ];
        [alphaest,cfCest,cfest,Cest,Q,Valpha,VcfC,CoValphacfC,Vcf,CoValphacf,unsafe,yj,varj,aest]= regrescomp(regu,nj,muj, j1,j2,0);
        alpha = [ alpha alphaest ];
        H  = [ H (alphaest + 1)/2 ];
        cf = [ cf cfest ];
        subQ = [ subQ Q ];
        stdalpha(n) = sqrt(Valpha) ;
        stdH(n) = stdalpha(n)/2;
        stdcf(n) = sqrt(Vcf);

        if (printout)
          if (wantH)
            fprintf('For subseries %d, (j1,j2)= (%d,%d), H = %4.2f, cf = %7.4f, Q=%5.4f .\n',n,j1,j2,H(n), cf(n),Q)
          else
            fprintf('For subseries %d, (j1,j2)= (%d,%d), alpha= %4.2f, cf = %7.4f, Q=%5.4f .\n',n,j1,j2,alpha(n),cf(n),Q)
          end 
          %%%  plot the next logscale diagram (in a row.)
          subplot(1,nsub,n)
          plot(yj,'*-')
          hold
          grid
          title(['Q = ',num2str(Q)])
          %title(['\alpha = ',num2str(alphaest,2),',  Q = ',num2str(Q)])
          %xlabel('Octave j')
          %ylabel('y(j)')
          upper = yj+seuil*sqrt(varj);
          lower = yj-seuil*sqrt(varj);
          jj = j1:j2;
          maxupperfull = max(maxupperfull,max(upper));
          minlowerfull = min(minlowerfull,min(lower));
          maxupper = max(maxupper,max(upper(jj)));
          minlower = min(minlower,min(lower(jj)));
          for k=jj    % plot vertical confidence intervals
             plot([k k ],[lower(k) upper(k)]) ;
          end
          plot(jj,alphaest * jj + aest,'r')
          %V = axis;
          %axis([0.8 j2+0.2  V(3) V(4) ])
          hold off
        end
      end
   end
   %  set LD's to the same vertical scale:
   if printout 
      for n = 1:nsub
        subplot(1,nsub,n)
        axis([j1-0.2 j2+0.2 minlower maxupper]);
        %axis([0.2 length(nj)+0.6 minlowerfull-0.5 maxupperfull-1.5]);% special for slide: full j range
      end
   end
   CIalpha = alpha;
   CIalpha(CIalpha<0)=0.01; %  if alpha -ve set to 0
   CIalpha(CIalpha>1)=.9;   %  if alpha >1  set to 1   
   cgam = cf./( 2*((2*pi).^CIalpha) .* gamma(CIalpha) .* sin(pi*(1-CIalpha)/2) );    %  need this for IC's for mean
   LRDasympstd = sqrt(cgam)./( lensub.^((1-CIalpha)/2) .* sqrt((1+CIalpha).*CIalpha/2) );  % asymptotic LRD CI's
   
   %%% Calculate and store H and cf over the entire series,  and the corresponding asympstd for the mean CI.
   %  Initialize the MRA based recursive method of calculating the DWT coefficients, then get the details 
   if  discrete_init==1    % use the special initialisation for intrinsically discrete series
      filterlength = 0;    %  choose the automatic length selection algorithm, or set here if desired
      [appro,kfirst,klast] = initDWT_discrete(x,regu,filterlength,0);      % no output
      if printout
         fprintf('** Using initialization on the full series, filterlength = %d\n',len+kfirst-klast)
      end
      %n = klast-kfirst+1;
      [muj,nj]=wtspec(appro,regu,fix( log2(klast-kfirst+1) ));
      clear appro                                       % this can be big! 
   else
      if printout
         fprintf('** Taking the given data as the initial approximation sequence\n')
      end
      [muj,nj]=wtspec(x,regu,fix( log2(len) ));
   end