function fw=FIEGARCH_spectral_density(lambda,parameters,ErrorType,zeta)
%%%%%%%%%%%%
%   USAGE %
%%%%%%%%%%%%
% Calculate the spectral density of a FIEGARCH model as a function of
% ordinary frequency, not angular frequency
% fw(j)/2pi=f(lambda);fw(j)=2pi*f(lambda);
% the autocovariance function can be then recovered by integrating the
% spectral density over interval <-1/2,1/2>

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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:FIEGARCH_spectral_density(lambda,parameters,ErrorType,zeta), [ErrorType,zeta] optional, stdnormal by default');
end

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

%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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     %
%%%%%%%%%%%%%%%%%%%%%%%%%%
% See Perez, Zaffaroni(2008)

if strcmp(ErrorType,'stdnormal')
                muz=sqrt(2/pi);
                alpha=psi(1,0.5);                    % psi(1,x)trigamma function
                gama=delta*muz*(psi(1)-psi(0.5));    % psi(x) 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);         
                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);     
                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);

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

cosl=cos(2*pi*lambda);         
omcosl=2*(1-cos(2*pi*lambda)); 

% Calculate Spectral Density Function fw 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

fw=(alpha+beta*fw1+gama*fw2); %/(2*pi):  transform from spectral density 
                              %defined for angular frequencies to transform 
                              %defined on simple frequencies (nx1)




%%%%%%%%%%%%%%%%
%   COPYRIGHT  %
%%%%%%%%%%%%%%%%
% Original code for Fourier-based Whittle: 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

%%%%%%%%%%%%%%%%%%%
%   REFERENCES    %
%%%%%%%%%%%%%%%%%%%
% PERCIVAL WALDEN: Wavelet Methods for Time Series Analysis
% Perez, Zaffaroni(2008): Finite-sample Properties of Maximum
%       Likelihood and Whittle Estimators in EGARCH and FIEGARCH Models