function  f = GPDEGARCHLogLik(x);
global GARCH_SERIES;
% GPD_EGARCH_LOG_LIK.M        Constructs the objective function involved 
%   in the ML estimation of the parameters of a GPD-EGARCH(1,1) process: 
%       Y[t] = mu + sig[t]*eta[t] where
%       log(sig[t]^2) = alpha0 + alpha1*log(sig[t-1]^2) +
%       alpha2*(Y[t-1]-mu) + alpha3*abs(Y[t-1]-mu)
%   and eta[t] is IID GPD distributed with mean zero, unit variance and pdf as given 
%   in Equation (3) in Komunjer (2006) with alpha = 0.5:
%       p(x,lda) = [gamma(2/lda)/(gamma(1/lda)*gamma(1+1/lda))]*[1/mu]*...
%                   exp(-[gamma(2/lda)/(alpha*gamma(1/lda))]^lda*[abs(x/mu+1-2*alpha)]^lda), if x <= -1*(1-2*alpha)*mu
%                   exp(-[gamma(2/lda)/((1-alpha)*gamma(1/lda))]^lda*[abs(x/mu+1-2*alpha)]^lda), if x > -1*(1-2*alpha)*mu
%       where mu = gamma(2/lda)*[gamma(3/lda)*gamma(1/lda)*(1-3*alpha+3*alpha^2)-gamma(2/lda)^2*(1-2*alpha)^2]^(-0.5).
%
% References:
%   [1] Komunjer, I. (2006): "Asymmetric Power Distribution: Theory and 
%   Applications to Risk Measurement", Journal of Applied Econometrics (forthcoming)
%
% Author: Ivana Komunjer - UCSD (komunjer@ucsd.edu)
% Version: 2.1        Date: 09/11/2006

alpha0 = x(1,1); alpha1 = x(2,1); alpha2 = x(3,1);  alpha3 = x(4,1); muY = x(5,1);% EGARCH parameters
a = 0.5; lda = x(7,1);  % alpha-quantile parameters
  
% Constructing the Innovations
Yt = GARCH_SERIES;
epst = Yt - muY;
T = size(Yt,1);
sig2t = (std(Yt)^2)*(1:1:T)';
for i = 2:1:T
    sig2t(i) = exp(alpha0 + alpha1*log(sig2t(i-1)) + alpha2*(epst(i-1)) + alpha3*abs(epst(i-1)));
end   
Innovations = epst./(sig2t.^0.5);
% we discard first 5%T values 
Innovations = Innovations(floor(0.05*T):T);
sig2t = sig2t(floor(0.05*T):T);
Tinn = size(Innovations,1);

% Computing the (minus) GPD log-likelihood for FMINCON
mu = gamma(3/lda)*gamma(1/lda)*(1-3*a+3*a^2)-((gamma(2/lda))^2)*(1-2*a)^2;
mu = gamma(2/lda)*(mu^(-0.5));
S = ((abs(Innovations/mu+1-2*a))*((gamma(2/lda))/(gamma(1/lda))));
S = (S.^lda).*[(Innovations<=-1*(1-2*a)*mu)/(a^lda) + (Innovations>-1*(1-2*a)*mu)/((1-a)^lda)];
S = sum(S);
f = -1*log((gamma(2/lda))/(gamma(1/lda)*(gamma(1+1/lda))*mu)) + 0.5*sum(log(sig2t))/Tinn + S/Tinn;

