%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   WAVELET WHITTLE: FIEGARCH(1,d,2)     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function WWLL=FIEGARCH_L_WW_estimate(method,wavelet,level,parameters,Ydata,ErrorType,zeta);

%%%%%%%%%%%%
%   USAGE %
%%%%%%%%%%%%
% Calculate the wavelet based (time-frequency) counterpart (WWLL) of the
% negative log-likelihood function (LL) with respect to restricted
% parameters of an FIEGARCH model and various error distributions. Various
% types of wavelet decomposition available, dealing with jumps in the time 
% series under analysis available.

%%%%%%%%%%%%%%%%
%   COPYRIGHT  %
%%%%%%%%%%%%%%%%
%   Lucie Kraicova
%   kraicova.l@seznam.cz

%%%%%%%%%%%%%%%%%%%
%   REFERENCES    %
%%%%%%%%%%%%%%%%%%%
% PERCIVAL WALDEN: Wavelet Methods for Time Series Analysis

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   INPUT CHECK AND DEFAULT VALUES FOR OPTIONAL PARAMETERS %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   Entries options, zeta decomposition :
%   Not used in montecarlo for master thesis


switch nargin
    case 5
        ErrorType='stdnormal';
        %zeta=[];
        nu=2;
        disp('Zt N(0,1) distributed')
    case 7
         if strcmp(ErrorType,'stdnormal')
                %zeta=[];
                nu=2;
                
         elseif strcmp(ErrorType,'stdged')
                 nu=zeta(1);
                 if nu<1
                     error('invalid nu')
                 end
                 
         elseif strcmp(ErrorType,'stdstudent')
                 nu=zeta(1);
                 if nu<=2
                     error('invalid nu')
                 end
                 
         elseif strcmp(ErrorType,'skewstudent')
                 nu=zeta(1);
                 if nu<=2 
                     error('invalid nu')
                 end
                 skew=zeta(2);
                 if abs(skew)>0.99
                     error('invalid lambda')
                 end
                 
         else
             error('invalid ErrorType')
         end
    otherwise
        error('5 or 7 inputs only:WWLL=FIEGARCH_L_W(wavelet,level,parameters,Ydata,ErrorType,zeta), [ErrorType,zeta] optional, stdnormal by default');
end

%}

%disp(['parameters=', num2str(parameters)]);

%%%%%%%%%%%%%%%%%%%%%%%%%%
%   INPUT:  DATA         %
%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Ydata= data assumed to be generated by FIEGARCH(q,d,p)*Nx1

T=length(Ydata);

%%%%%%%%%%%%%%%%%%%%
%  INPUT: METHOD   %
%%%%%%%%%%%%%%%%%%%%
% Methods Available:
%   'DWT' (uses wmtsa toolbox)
%       Fixed parameters:
%           boundary   -- (optional) boundary conditions to use (string)
%                   Valid values: 'circular' or 'reflection'
%                   Default: 'reflection'
boundary='circular';
%           * opts       -- (optional) Additional function options
%              Optional input arguments are specified as name-value pairs:
%                   * RetainVJ -- Boolean flag indicating whether to
%                           scaling coefficients have been retained at all
%                           levels.
%                       Values:  1 = true,  all VJ retained,
%                                0 = false, VJ retained for J0 level.
%                       Default: 0, VJ retained only at J0 level.
opts='';
        

%   'dwt' (uses matlab toolbox (function wavedec))


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                      Nothing else to be specified                      %
%(only inclusion of boundary coefs. and of scaling coefs - not default   %
% and setting appropriate methods for desling with jumps)                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%
%  INPUT: WAVELETS %
%%%%%%%%%%%%%%%%%%%%
% Filters Available:
%   For WMTSA transform(MODWT): write: 
%       wtfilter('list')
%       or a valid filter struct

%   For MATLAB transform(DWT,SWT,NDWT): write: 
%       waveletfamilies('n')
%       or a valid filter struct

%disp(['wavelet= ', num2str(wavelet)]);

%%%%%%%%%%%%%%%%%%%%
%   INPUT: LEVEL   %
%%%%%%%%%%%%%%%%%%%%
%   0 < integer < log2(T)
%   'conservative'  <log2( (T / (L - 1)) - 1) (L is the wavelet support)
%   'max' <log2(T)
%   'supermax' <log2(1.5 * T)

%  'REFLECTION'
%   'dwtcut3'  <floor(log2(T))-3  16 J-level coefficients in DWT
%   'dwtcut4'  <floor(log2(T))-4  32 J-level coefficients in DWT
%   'dwtcut5'  <floor(log2(T))-5  64 J-level coefficients in DWT
%   'dwtcut6'  <floor(log2(T))-6 128 J-level coefficients in DWT
%   'dwtcut7'  <floor(log2(T))-7 256 J-level coefficients in DWT
%   'dwtcut8'  <floor(log2(T))-8 512 J-level coefficients in DWT

%  'CIRCULAR'
%   'dwtcut3'  <floor(log2(T))-3  8 J-level coefficients in DWT
%   'dwtcut4'  <floor(log2(T))-4  16 J-level coefficients in DWT
%   'dwtcut5'  <floor(log2(T))-5  32 J-level coefficients in DWT
%   'dwtcut6'  <floor(log2(T))-6  64 J-level coefficients in DWT
%   'dwtcut7'  <floor(log2(T))-7 128 J-level coefficients in DWT
%   'dwtcut8'  <floor(log2(T))-8 256 J-level coefficients in DWT

% Translate the level input into valid integer

% Not useful in case of modwt with boundary coefficients restrictions

%{
    %for strings:
if strcmp(level,'max')
    level=floor(log2(T));
    %disp(['level set to maximum: ','level= ', num2str(level)]);
 elseif strcmp(level,'supermax')
    level=floor(log2(1.5 * T));
    %disp(['level set to supermax value: ','level= ', num2str(level)])
 elseif strcmp(level,'conservative')
    if strcmp(method,'modwt')||strcmp(method,'DWT')
        filter = modwt_filter(wavelet);
        L=filter.L;
    elseif strcmp(method,'dwt')||strcmp(method,'swt')||strcmp(method,'ndwt')
        [LB, UB] = wavsupport(wavelet);
        L=UB-LB+1;
    else
        error('method not known to level translator: add it or choose different level')
    end
    level=floor(log2( (T / (L - 1)) - 1));
    %disp(['level set to conservative value: ','level= ', num2str(level)])
 elseif strcmp(level,'dwtcut3')
    level=floor(log2(T))-3;
    if level<1
        level=floor(log2(T));
        disp(['not enough observations; level reduced to maximal possible value: ','level= ', num2str(level)])
    else
        %disp(['level set to dwtcut value: ','level= ', num2str(level)]);
    end
 elseif strcmp(level,'dwtcut4')
    level=floor(log2(T))-4;
    if level<1
        level=floor(log2(T));
        disp(['not enough observations; level reduced to maximal possible value: ','level= ', num2str(level)])
    else
        %disp(['level set to dwtcut value: ','level= ', num2str(level)]);
    end
 elseif strcmp(level,'dwtcut5')
    level=floor(log2(T))-5;
    if level<1
        level=floor(log2(T));
        disp(['not enough observations; level reduced to maximal possible value: ','level= ', num2str(level)])
    else
        %disp(['level set to dwtcut value: ','level= ', num2str(level)]);
    end
 elseif strcmp(level,'dwtcut6')
    level=floor(log2(T))-6;
    if level<1
        level=floor(log2(T));
        disp(['not enough observations; level reduced to maximal possible value: ','level= ', num2str(level)])
    else
        %disp(['level set to dwtcut value: ','level= ', num2str(level)]);
    end
 elseif strcmp(level,'dwtcut7')
    level=floor(log2(T))-7;
    if level<1
        level=floor(log2(T));
        disp(['not enough observations; level reduced to maximal possible value: ','level= ', num2str(level)])
    else
        %disp(['level set to dwtcut value: ','level= ', num2str(level)]);
    end     
 elseif strcmp(level,'dwtcut8')
    level=floor(log2(T))-8;
    if level<1
        level=floor(log2(T));
        disp(['not enough observations; level reduced to maximal possible value: ','level= ', num2str(level)])
    else
        %disp(['level set to dwtcut value: ','level= ', num2str(level)]);
    end      
  %for integers:
 elseif level>floor(log2(T))
    level=floor(log2(T));
    disp(['level reduced to maximal possible value: ','level= ', num2str(level)])
 elseif level>0 && level<=floor(log2(T))
    %disp(['level is valid: ','level= ', num2str(level)])
 else
    level=floor(log2(T));
    disp(['invalid level input -> level set to maximal possible value: ','level= ', num2str(level)])
 end

%} 

%%%%%%%%%%%%%%%%%%%%%%%%%%
%   INPUT:  PARAMETERS   %
%%%%%%%%%%%%%%%%%%%%%%%%%%
%parameters=[d,omega,aL,bL,theta,delta]

d=parameters(1);
omega=parameters(2); %not available for both Fourier-bsed and wavelet whittle 
a2=parameters(3);          
b1=parameters(4);                
theta=parameters(5);        
delta=parameters(6); 

theta_2=theta^2;


% Restrictions:
    %  -0.5<d<0.5
    %  omega=0 assumption 
	%  -1<a1<1  
    %  -1<b1<1             
    % 0<theta_2
    % delta finite

%   FREQUENCIES:

Frequencies=1./(2.^((1:1:level+1)));
Frequencies = Frequencies';
% Respective angular frequencies= 2pi*Frequencies

delta_DIRAC=ones(level+1,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       WAVELET TRANSFORM       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%{
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                     JUMP DETECTION                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
[WJt, VJt, att] = modwt(Ydata,wavelet,1,boundary,opts);  % Modwt transform
[WJt, VJt] = modwt_cir_shift(WJt, VJt, wavelet, 1); % Circular shift to match true indices
Detail=WJt(1:size(WJt),1); % Get wavelet coefficients at first level
           
  % Simple Jump location detection: Hard thresholding on first level coeffs
        % Donoho and Johnstone 1994 threshold
         thrs=(2*log(length(Detail)))^(0.5)*2^(0.5)*median(abs(Detail))/0.6745;                                                     
          
         jump_loc=abs(Detail)>thrs;
         nojump_loc=abs(jump_loc-1);
               
         Yavg=[(Ydata(end)+Ydata(2))/2;(Ydata(1:end-2)+Ydata(3:end))/2;(Ydata(end-1)+Ydata(1))/2]; % jump correction for all locations
         
         Ydata=nojump_loc.*Ydata+jump_loc.*Yavg; % do correction at jump_loc indices         
 
 %}        
         
         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                     WWE                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xt=log(Ydata.^2);   
xtc=xt; %-mean(xt); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                     CHOOSE TYPE OF WAVELET TRANSFORM                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% [1]   USING MATLAB CODE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Coef.Variance: dwt     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if strcmp(method,'dwt')
   % not adjusted to exclude boundary coefficients
   % display('Type of wavelet transform= DWT (Discrete Wavelet Transform)')
    [WavCoefs,CoefsLength]=wavedec(xtc,level,wavelet); 
    Details=detcoef(WavCoefs,CoefsLength,[1:level]);

 %Get coefficients at individual levels and compute sums of their squares
    coefcount_vec=zeros(level+1,1);%Vector of numbers of coefficients at each level
    Coef2_vec=zeros(level+1,1);
    for j=1:level
        Detail=cell2mat(Details(j)); % No substraction of boundary coefficients here
        coefcount_vec(j)=length(Detail);
        Coef2_vec(j)=Detail'*Detail;
    end
    
    % Using scaling coefficients
    
    %{
    scaling_coef=appcoef(WavCoefs,CoefsLength,wavelet);
    coefcount_vec(level+1)=length(scaling_coef);
    %}
    % For processes with non-zero mean
    %{
     scaling_coef=scaling_coef-mean(scaling_coef); 
    %}
    
    %{
    Coef2_vec(level+1)=scaling_coef'*scaling_coef;
    %}
    

    
% [2]   USING PERCIVAL WALDEN CODE FOR "Wavelet Methods for Time Series Analysis
%       Textbook"   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Coef.Variance: DWT     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(method,'DWT')
    %display('Type of wavelet transform= DWT (Discrete Wavelet Transform)')
    
    [~, ~, LJ] = dwt_equivalent_filter(wavelet, level); % lengths of rescaled filter LJ if boundary coefficients should be ommited (first Lj/2)
    [WJ, VJ, ~, NJ] = DWT(xtc,wavelet,level,boundary,opts);

%Get coefficients at individual levels and compute sums of their squares
   
    Coef2_vec=zeros(level+1,1);
    coefcount_vec=zeros(level+1,1);
    for j=1:level
        Detail=cell2mat(WJ(j));
        Detail=Detail(cell2mat(LJ(1))-1:length(Detail)); %max LJ(1)-2 boundary coeffs. (based on Percival Walden 2000)
        coefcount_vec(j)=length(Detail);
        
        Coef2_vec(j)=Detail'*Detail;   
    end
    % Using scaling coeffs
    %{
    scaling_coef=cell2mat(VJ(length(VJ)));
    scaling_coef=scaling_coef(cell2mat(LJ(1))-1:length(scaling_coef));
    coefcount_vec(level+1)=length(scaling_coef);
    %}
    
    %For processes with non-zero mean  
    %{
    %scaling_coef=scaling_coef-mean(scaling_coef);
    %}
    
    %{
    Coef2_vec(level+1)=scaling_coef'*scaling_coef;
    %}
    
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Coef.Variance: MODWT     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif strcmp(method,'modwt')
    %display('Type of wavelet transform= MODWT (Maximal Overlap Discrete Wavelet Transform)')
    
    [hJt, ~, LJ] = modwt_equivalent_filter(wavelet, level); % lengths of rescaled filter LJ if boundary coefficients should be ommited (first Lj-2)
    [WJt, VJt, att] = modwt(xtc,wavelet,level,boundary,opts); 
    

%Get coefficients at individual levels and compute sums of their squares
    
    Coef2_vec=zeros(level+1,1);
    coefcount_vec=zeros(level+1,1);%Vector of numbers of coefficients at each level
    coefcount_restr=zeros(level+1,1);%Vector of numbers of coefficients at each level
    
    for j=1:level
        Detail=WJt(1:size(WJt),j);  
        coefcount_vec(j)=length(Detail); % number of coefs equal to the sample length
        Detail=Detail(cell2mat(LJ(j)):length(Detail));% First LJ(j)-1 coeffs 
             % are boundary coeffs (Should be based on Percival and Walden
             % 2000)        
        coefcount_restr(j)=length(Detail);%Vector of numbers of coefficients 
                      % at each level after substracting boundary coeffs used 
                      % for wavelet variance estimation variance                    
         
         Coef2_vec(j)=Detail'*Detail;
    end

    
       % Using scaling coeffs
    %{
    scaling_coef=VJt;    
    coefcount_vec(level+1)=length(scaling_coef);
    scaling_coef=scaling_coef(cell2mat(LJ(j))-1:length(scaling_coef)); 
    coefcount_restr(level+1)=length(scaling_coef); 
   
    %For processes with non-zero mean    
      %{
       scaling_coef=scaling_coef-mean(scaling_coef); 
      %}
    
    Coef2_vec(level+1)=scaling_coef'*scaling_coef;
    %}
    
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  OBJECTIVE FUNCTION       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

spectral_density=@(lambda)FIEGARCH_spectral_density(lambda,parameters,ErrorType,zeta);

fw=zeros(level+1,1);
  for j=1:level
    fw(j)=integral(spectral_density,Frequencies(j+1),Frequencies(j),'AbsTol',1e-7 );
  end
    %fw(level+1)=integral(spectral_density,1/(500000000),Frequencies(level+1));

  for k=1:level+1
    if isnan(fw(k)) || isinf(fw(k))
    WWLL=1.0e+20;
    NAN=1;
    else
    NAN=0;
    end
  end

 if  NAN==0;
   % Using both wavelet and scaling coefficients (enable calculation of
      % fw(level+1) for data with non-zero mean enable its substraction in
       % Coef2_vec calculation )
   
      %Frequencies(level+1)=Frequencies(level)
      %WWLL=(coefcount_vec(1:level+1))'*log(2.*fw(1:level+1)./Frequencies(1:level+1))+(delta_DIRAC(1:level+1))'*(Coef2_vec(1:level+1).*Frequencies(1:level+1)./(2.*fw(1:level+1)));    
  
   
   if strcmp(method,'dwt') || strcmp(method,'DWT')
     WWLL=(coefcount_vec(1:level))'*log(2.*fw(1:level)./Frequencies(1:level))+(delta_DIRAC(1:level))'*(Coef2_vec(1:level).*Frequencies(1:level)./(2.*fw(1:level))); 
   
   elseif strcmp(method,'modwt')
     WWLL=(coefcount_vec(1:level).*Frequencies(1:level))'*log(2.*fw(1:level)./Frequencies(1:level))+(delta_DIRAC(1:level))'*(coefcount_vec(1:level)./coefcount_restr(1:level).*Coef2_vec(1:level).*Frequencies(1:level)./(2.*fw(1:level)));   
   
   end
   
 end

end



