%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   MAXIMUM LIKELIHOOD: FIEGARCH(1,d,2)  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [LL2estimate]=FIEGARCH_L_ML_estimate(parameters,Ydata,ErrorType,zeta)

%disp('[LL2estimate]=FIEGARCH_L_ML_estimate(parameters,Ydata,ErrorType,zeta)')
%disp(['parameters=', num2str(parameters)]);

%%%%%%%%%%%%
%   USAGE %
%%%%%%%%%%%%
% Returns value of negative log-likelihood and keeps estimation working
% even in case of anomalous situations.

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

%   Entries options, zeta decomposition :

switch nargin
    case 2
        ErrorType='stdnormal';
        zeta=[];
        %disp('Zt N(0,1) distributed')
    case 4
         if strcmp(ErrorType,'stdnormal')
                zeta=[]; 
                 
         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 %disp('invalid ErrorType')
         end
    otherwise
        error('2 or 4 inputs only:[Y,ht,Phi]=FIEGARCH_simulate(M,N,d,omega,aL,bL,theta,delta,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

%%%%%%%%%%%%%%%%%%%%%%%%%%
%   INPUT:  PARAMETERS   %
%%%%%%%%%%%%%%%%%%%%%%%%%%
%parameters=[d,omega,aL,bL,theta,delta]
d=parameters(1);
omega=parameters(2);
a2=parameters(3);          
b1=parameters(4);                
theta=parameters(5);        
delta=parameters(6);  

%Restrictions:
    %  -0.5<d<0.5
    %  omega finite
	%  -1<a2<1  
    %  -1<b1<1             
    %  theta finite 
    %  delta finite 
    
    

%{
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                     JUMP DETECTION                                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wavelet='d4'; % Filters Available:wtfilter('list') 
%be careful about the number of coefficients the jump influences (Haar:2-symetrically, D4:3 - sum of 2 smaller effects=-bigger effect)
boundary='circular'; 
opts='';

[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  
       
%}
   
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  NEGATIVE LOG-LIKELIHOOD %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Calculate vector of conditional variances
[Yt,eht,Zt]=FIEGARCH_core(Ydata,d,omega,a2,b1,theta,delta, ErrorType,zeta);

%Do not call distr.m in cases of invalid conditional variances, set extreme
%LL value instead-> extremely high neg.log-likelihood, it keep the
%estimation going

if isnan(eht(length(eht)))  
    LL=-1.0e+20;
elseif eht(length(eht))==0  
    LL=-1.0e+20;
elseif isinf(eht(length(eht)))  
    LL=-1.0e+20;
else
    % avg=mean(Yt);
     %Calculate negative Log-likelihood functions for any valid distribution:
    if strcmp(ErrorType,'stdnormal');  
      % disp('N(0,1) errors');
        LL= normloglik(Yt,0,eht);
        
    elseif strcmp(ErrorType,'stdged')  
       % disp('Standardized GED(nu) errors');
        LL = gedloglik(Yt,0,eht,nu);
           
     elseif strcmp(ErrorType,'stdstudent')  
     % disp('Standardized Student-t(nu)errors');
        LL = stdtloglik(Yt,0,eht,nu);
         
     elseif strcmp(ErrorType,'skewstudent') 
     % disp('Standardized Skewed Student-t(nu, lambda) errors');
        LL= skewtloglik(Yt,0,eht,nu,skew);
        
     end;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Should not be important, see correction below FIEGARCH_core

% Catch conditions that produce anomalous log-likelihood function values.
% Typically, what happens is that input parameter values will result in
% an unstable inverse filter, which produces LLF = inf. This, in turn, will
% not allow FMINCON to update the iteration. By setting the LLF to a large, 
% but finite, value, we can safeguard against anomalies and allow the 
% iteration to continue. (correction used by Ana Perez)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%{
if isinf(LL)
    LL  =  1.0e+20;
elseif isinf(-LL)
    LL=-1.0e+20;
elseif isnan(LL)
    LL=1.0e+20;
end
%}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%disp(['LL=',num2str(LL)])
LL2estimate=-LL;

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

% Toolboxes used: MFE Toolbox by Kevin Sheppard
%   http://www.kevinsheppard.com/wiki/MFE_Toolbox; downloaded 19/02/2013: 
%           gegloglik.m (modified to allow nu=1)
