%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   WHITTLE: FIEGARCH(1,d,2)             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function WLL=FIEGARCH_L_W_estimate(parameters,Ydata,ErrorType,zeta);

%%%%%%%%%%%%
%   USAGE %
%%%%%%%%%%%%
% Calculate the frequency domain counterpart (WLL) of the negative
% log-likelihood function (LL).

%%%%%%%%%%%%%%%%
%   COPYRIGHT  %
%%%%%%%%%%%%%%%%
% Original code: Ana Perez, code for paper: Perez,
%   Zaffaroni(2008),Finite-sample Properties of Maximum Likelihood and
%   Whittle Estimators in EGARCH and FIEGARCH Models, received on request
%   25/01/2013 

% Modifications: Lucie Kraicova (kraicova.l@seznam.cz) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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
             error('invalid ErrorType')
         end
    otherwise
        error('2 or 4 inputs only:[WLL,WLLS]=FIEGARCH_L_W(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

%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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); 

theta_2=theta^2;

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

%%%%%%%%%%%%%%%%%%%%%%%%%%
%   SPECTRAL DENSITY     %
%%%%%%%%%%%%%%%%%%%%%%%%%%

 if strcmp(ErrorType,'stdnormal')
                muz=sqrt(2/pi);
                alpha=psi(1,0.5);                    % trigamma function
                gama=delta*muz*(psi(1)-psi(0.5));    % digamma function
                
 elseif strcmp(ErrorType,'stdged')
                muz=gamma(2/nu)/sqrt(gamma(3/nu)*gamma(1/nu));
                alpha=(2/nu)^2*psi(1,1/nu);          % trigamma function
                gama=2*(delta/nu)*muz*(psi(2/nu)-psi(1/nu));
                
                 
 elseif strcmp(ErrorType,'stdstudent')
                muz=sqrt(nu-2)*gamma((nu-1)/2)/(sqrt(pi)*gamma(nu/2));
                alpha=psi(1,0.5)+psi(1,nu/2);      % trigamma function
                gama=delta*muz*(psi(1)-psi((nu-1)/2)-psi(0.5)+psi(nu/2));
               
                 
 elseif strcmp(ErrorType,'skewstudent')
          error('L_W and L_WW not specified for this distribution')
                 
 else disp('invalid ErrorType') 
 end

beta=theta_2+delta^2*(1-muz^2);


T=length(Ydata);
n=fix(T/2);                     % integer part of T/2 for periodogram
delta_DIRAC=ones(n,1);
 if n == (T/2);                  % cheking whether T is even 
    delta_DIRAC(n)=1/2;
 end;

lambda=(2*pi*(1:1:n)')/T; % nx1 because of the periodogram

S1=real((1-exp(i*lambda)).^(-d) );
S2=real( exp(i*lambda).*(1-exp(i*lambda)).^(-d) );
S3=real( exp(i*2*lambda).*(1-exp(i*lambda)).^(-d) );

cosl=cos(lambda);         
omcosl=2*(1-cos(lambda)); 
%disp(['max[S1trun-S1,S2trun-S2,S3trun-S3] = ', num2str(max([S1trun-S1 S2trun-S2 S3trun-S3]))]);

% Calculate Spectral Density Function of Xt from FIEGARCH(1,d,2)
fw1=(omcosl.^(-d)).*(1+2*a2*cosl+a2^2)./(1-2*b1*cosl+b1^2); % nx1
fw2=2*((1-a2*b1)*S2+a2*S3-b1*S1)./(1-2*b1*cosl+b1^2);       % nx1

% Spectral density = fw =(alpha+beta*fw1+gama*fw2)/(2*pi);   nx1 vector
% the autocovariance function can be recovered by integrating the
% spectral density over interval <-pi,pi>

% In the objective function: 
    % 1] log(fw)=log(fw.*1)=log(fw.*2pi/2pi)=log(fw.*2pi)-log(2pi) the later
        % term can be ommited in the objective function, its a constant
    
    % 2] Periodogram(resp.to 2pij/T)/fw(2pij/T):
        % fw(j), j=2,..N,respective frequencies f=1/N,...(N-1)/N, angular frequencies lambda=2pi*1/N,...2pi*(N-1)/N
        % fw(j)/2pi=f(lambda);f(j)=2pi*f(lambda); thus 2pi in the spectral
        % density cancels out with the 2pi which rescale the periodogram.
        % To make computation somewhat faster, together with [1] I can use
        % adjusted spectral density and adjusted periodogram.
        
        
fw=(alpha+beta*fw1+gama*fw2);%/(2*pi);




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

 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                     FWE                                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%
%      PERIODOGRAM       %
%%%%%%%%%%%%%%%%%%%%%%%%%%

xt=log(Ydata.^2);   
xtc=xt-mean(xt);   
periodog=(abs(fft(xtc))).^2/(length(xtc)); %/2*pi (adjusted as discussed above) % Tx1 (I0,I1,...,IT-1)    
I_vec=periodog(2:n+1);            % nx1 column vector with I1,...,I[T/2] 
clear xtc periodog;

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

WLL=delta_DIRAC'*log(fw)+delta_DIRAC'*(I_vec./fw);   
end

