%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                                                       %%   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 variableswantH = 0;   % print and output according to H or alpha,   alpha = 2H - 1%%%%%%%%%%% correct inputif  (lookH | lookmean | lookcf)  % if want to look at H, make sure regu is set.   regu=max(regu,1);enddim = size(x);if dim(1) > 1  x = x';end%%%%%%%%%%% initx = x';len = length(x);m = 0;v = 0;al = 0;stdal = 0;fullal = 0;fullstdal = 0;crit_level = 0;%%%%%%%%%%%  determine subseries of xif (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 A4else  nsub = length(cutpoints) + 1;  lensub = diff([1 cutpoints len]) % *** after this unequal lengths are not supported! but could be doneend    %%%%%%%%%%  create subseries as rows of matrix 'series'series = ones(nsub,lensub);for n = 1:nsub  series(n,:) = x( (1+(n-1)*lensub):(n*lensub) )';endseries(:,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 subseriesseuil = 1.9599 ;num_plot = ( looksig + lookmean + lookvar + lookH + lookcf );nplotnum = (num_plot)*100 + 10;if (printout)  figure(20)  clfendif (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')endline   = 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