%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                                                       %%   LDestimate.m                                        %%                                                       %%        D. Veitch   P.Abry                             %%                                                       %%   1/6/98                                              %%  DV 4/99                                              %%  DV 7/99                                              %%  DV 15/5/2000                                         %%  DV 8/10/2002 Add extra outputs yj and varj           %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    This function estimates the scaling parameter  'alpha'  of a scaling process.%    Such processes include  Long Range Dependent (LRD), Self-Similar, 1/f noise, fractal,%    and multifractal processes.   In the case of LRD a second parameter, cf,  is also estimated,%    as LRD requires (at least) two parameters for its description.%      Note:  This function only looks at second order statistics (essentially covariances), this will not give%             the full picture for multifractal processes.%%    More precisely, this function estimates the two parameters of LRD: (alpha,cf), using the%    wavelet based joint estimator of Abry and Veitch.  It is a spectral domain definition, namely where%    the spectrum takes the form%                                          f(nu) ~ cf*nu^alpha near the origin, for a certain range of%    scales (frequencies).%%    In the case where the series to be analysed is intrinsically discrete, a special initialisation procedure%    must be followed to avoid errors at low scales.  This is activated via one of the arguments.%% *** Usage:  [alphaest,cfCest,cfest,Cest,Q,j1opt,yj,varj] = LDestimate(data,regu,j1,j2,discrete_init,calcj1,printout)%%--- Routines called Directly:% wtspec.m%      " function [muj,nj]=wtspec(data,regu,nbvoies) "% regrescomp.m%      " function [alphaest,cfCest,cfest,Cest,Q,Valpha,VcfC,CoValphacfC,Vcf,CoValphacf,unsafe]%                                                =  regrescomp(regu, nj, muj,j1,j2,printout) "% initDWT_discrete.m%      " function [appro,kfirst,klast] = initDWT_discrete(data,regu,lengthIwant,printout) "% gauss_CDFinv.m%      " function [quantile] = gauss_CDFinv(prob, mean, std) "% newchoosej1.m%     "  function [j1opt,Qmat] = newchoosej1(regu, nj, muj, printout, j2vec)"%%%  Input:   data:  the input data as a row vector:  this is in fact the sequence of wavelet "approximation"%                  coefficients, or an approximation thereof (for example often sampled data is used here).%           regu:  (regularity) =  number of vanishing moments (of the Daubechies wavelet).%           j1:    the lower limit of the scales chosen,  1<= j1 <= scalemax-1%           j2:    the upper limit of the octaves chosen, 2<= j2 <= scalemax%           discrete_init:   1:  perform the special MRA initialisation for intrinsically discrete series%                         else: assume input data already initialised (ie is already the approximation sequence)%           calcj1:  1: decision to run the newchoosej1 function, which will output a plot of Q(j_1) %                    vs j_1, as well as returning the optimal j1 values (one per j2 input to it).%                    The values of j2 chosen are not input parameters- need to edit this function.%                    NOTE: choosing this does Not affect the j1 passed to%                    LDestimate for estimation. %                 else: don't run newchoosej1 %           printout:   1:%                        -- a log-log plot graph is plotted plus the regression line and +- 1.96*std(muj)%                           around each point, being the 95% confidence interval under Gaussian assumptions%                        -- the values of H etc are printed.%                        -- a loop is entered allowing interactive choosing of  (j1,j2) (return to exit)%                    else:   nothing is printed or plotted.%%  Output:  alphaest:   estimate of the spectral estimate  alpha%           cfCest:     estimate of the intermediate quantity  cfC%           cfest:      estimate of the second LRD parameter  cf = cfC / C%           Cest:       estimate of the wavelet dependent integral C(alpha) using the estimated alphaest%           Q:          the goodness of fit measure, the Chi2 based probability that the null hyp is true given%                       the observed data over the scale range  (j1,j2)%           j1opt:      the optimal j1^*(j2) vector outputted from newchoosej1  %           yj:         the LD values%           varj:       their variances, allowing CI's to be replotted.%%   eg of interactive use%    >> load fgn8.dat                   % discrete data 4096 long, which is  2^12   (12 octaves)%    >> LDestimate(fgn8,3,1,12,1,1,1);    % try the full range initially, (j1,j2)=(1,12),  print and plot output%                                       % and initialise as the data is intrinsically discrete%                                       %  vary (j1,j2) interactively.%                                       % also run newchoosej1 to get its opinion%%   eg of batch use, with (j1,j2)=(4,11),  regu = 2  established by prior experiments.%   data is sampled real data, so do not performed the discrete initialisation%    for i= 1:no_realisations%     [alphaest[i],cfest[i] = LDestimate(realdata(i,:)),2,4,11,0,0,0);     % don't printout,%                                                     % only store (alpha,cf)%                                                     % don't initialize,%                                        % don't calculate optimal values as don't store them%    end%--------------------------------------------------------------------------------------------------------function [alphaest,cfCest,cfest,Cest,Q,j1opt,yj,varj] = LDestimate(data,regu,j1,j2,discrete_init,calcj1,printout)%%%  Initializeformat compact                          % eliminate excess blank lines in outputn = length(data);    % record original length%  Massage the input parametersnbvoies = fix( log2(n) );        % determine the largest possible number of octaves in the data (never realised)%  Initialize the MRA based recursive method of calculating the DWT coefficientsif  discrete_init==1    % use the special initialisation for intrinsically discrete series   filterlength = 0;    %  choose the automatic length selection algorithm, or set here if desired   %  here would be nice to call the output 'appro',  but wastes memory   [data,kfirst,klast] = initDWT_discrete(data,regu,filterlength,0);          % no output   if printout      fprintf('** Using initialization for discrete series, filterlength = %d\n',n+kfirst-klast)   end   n = klast-kfirst+1;else   if printout      fprintf('** Taking the given data as the initial approximation sequence\n')   endend%  perform the decomposition using Daubechies wavelets[muj,nj]=wtspec(data,regu,nbvoies);if printout  fprintf('No of points n_j at octave j:   ')  fprintf('%d ',nj)  fprintf('\n')  fprintf('Number predicted by nj=n*2^j:   ')  %fprintf('%d ',floor(nj(1).*2*2.^(-(1:length(nj)))) )  fprintf('%d ',floor(n*2.^(-(1:length(nj)))) )  fprintf('\n\n')end % Make an automated choice of j1  (value not used in the LD, but separate graph and text displayed)if (calcj1)   maxj2 = length(nj);   %[j1opt] = newchoosej1(regu, nj, muj, printout,4:maxj2);          % full range of possible j2   %[j1opt] = newchoosej1(regu, nj, muj, printout,[5,8 maxj2]);      % specific choice, optimise this   [j1opt] = newchoosej1(regu, nj, muj, printout, maxj2);            % LRD choice   %%% modify plot output of regresscomp   if discrete_init & printout      titlehandle = get(gca,'title');      titlestring = get(titlehandle,'string');   % get current title string      set(titlehandle,'string',[ titlestring, ',   D-init' ])  % append discrete init indication   endelse   j1opt = 1;   % return dummy value to avoid output errorend% **************** Interactive loop while  1 %  Verify input parameters j2 = min(j2,length(muj));    % make sure j2 cannot exceed the maximum number of octaves available j2 = max(j2,2);              % make sure that j2>1 j1 = min(j1,j2-1);           % make sure that j1<j2  %  Perform the joint parameter estimations and calculate the goodness of fit measure, and plot [alphaest,cfCest,cfest,Cest,Q,Valpha,VcfC,CoValphacfC,Vcf,CoValphacf,unsafe,yj,varj] = regrescomp(regu, nj,muj, j1,j2, printout);  %  Calculate the safe octaves for the cfC calculation, and printout a warning if printout=FALSE and some not safe  safe_octaves = j1:(j2-unsafe);  if  length(safe_octaves) ~= (j2-j1+1) & ~printout    fprintf('Warning, %d   octaves  were unsafe, not used to calculate cfC. \n',unsafe) end  % Printout a summary of the octaves: in-data/used/available/safe,   and the results on H etc if (printout)  %%% modify plot output of regresscomp   if discrete_init      titlehandle = get(gca,'title');       titlestring = get(titlehandle,'string');   % get current title string       set(titlehandle,'string',[ titlestring, ',   D-init' ])  % append discrete init indication   end   %divergence_test   %%% Text output:    fprintf('Octaves:    in_data     available     selected     Goodness of fit (Prob of data assuming lin-regression)\n')   fprintf('            1--%3.1f       1--%d          %d--%d                      %6.5f \n\n', log2(n), length(nj),j1,j2,Q )   %%%% Calculate Confidence Intervals,   95% two sided. Regresscomp outputs variance, here map to CI   %  set confidence level.     sig_level = 5 ;   % cf  95% two sided lognormal assumption.  L(M,V)=exp(N(m,v))   m = ln(M^4/(V+M*M))/2;  v = ln(V/M^2 +1)   z1 = sqrt(2) * erfinv(2*   sig_level/2/100  -1) ;   z2 = sqrt(2) * erfinv(2*(1-sig_level/2/100) -1) ;   m = log( (cfest^4)/(Vcf+cfest*cfest) )/2;   v = log( Vcf/cfest/cfest + 1 );   z1 = ( z1*sqrt(v) ) + m ;   z2 = ( z2*sqrt(v) ) + m ;    cfL = exp(z1);    cfR = exp(z2);    % Scaling parameters:  95% two sided gaussian assumption    seuil=1.9599;   HLRD = (alphaest + 1)/2;   H    = (alphaest - 1)/2;   h    = (alphaest - 1)/2;     D    = (5 - alphaest)/2;   aL = alphaest - seuil*sqrt(Valpha);   aR = alphaest + seuil*sqrt(Valpha);    HLRDL = HLRD - seuil*sqrt(Valpha/4);   HLRDR = HLRD + seuil*sqrt(Valpha/4);   HL    = H    - seuil*sqrt(Valpha/4);   HR    = H    + seuil*sqrt(Valpha/4);   DL    = D    - seuil*sqrt(Valpha/4);   DR    = D    + seuil*sqrt(Valpha/4);      %  Print the output   % fprintf('Goodness of fit statistic Q (probability of data assuming regression valid): %6.5f \n \n',Q)   fprintf('Scaling parameters are:    alpha (LRD)     H (LRD rewrite)  H=h (ss,Holder)  D (frac dim, if alpha in (1,3)) \n')   fprintf('             Estimates:       %4.3f            (%4.3f)           %4.3f            %4.3f\n',alphaest,HLRD,H,D)   fprintf('                  CI''s:  [%4.3f, %4.3f]   [%4.3f, %4.3f]   [%4.3f, %4.3f]   [%4.3f, %4.3f]\n\n', aL,aR,HLRDL,HLRDR,HL,HR,DL,DR )   fprintf(' Second parameters are:         cf              N/A            sigma^2              N/A         \n')   fprintf('             Estimates:      %7.4f                      Work in progress            \n',cfest)   fprintf('                  CI''s:  [%7.5f, %7.5f]                Work in progress          \n', cfL,cfR )   fprintf('\n\n')          %  Prompt for new (j1,j2) values   j1=input('New initial octave j1? (hit return to exit loop)   ');   if  isempty(j1)       hold off      return   end   j2=input('New final octave j2?   ');   hold off   fprintf('**************************************************************************************\n'); else   %  if you can't see the answer there's no point prompting for more values   return end hold offend