function  f = APDGARCHLogLik(x);
global GARCH_SERIES;
% APD_GARCH_LOG_LIK.M       Constructs the objective function involved in 
%   the ML estimation of the parameters of a GARCH(1,1) process 
%   with APD innovations:
%       Y[t] = mu + sig[t]*eta[t] where
%       sig[t]^2 = alpha0 + alpha1*sig[t-1]^2 + alpha2*(Y[t-1]-mu)^2 and
%       eta[t] is IID APD distributed with mean zero and unit variance 
% The pdf of eta[t] is as in Equation (3) in Komunjer (2006):
%       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

omega = x(1,1); alpha = x(2,1); beta = x(3,1);  muY = x(4,1);% GARCH parameters
a = x(5,1); lda = x(6,1);  % alpha-quantile parameters
  
% Constructing the Innovations
Yt = GARCH_SERIES;
epst = Yt - muY;
T = size(Yt,1);
sig2t = omega/(1-alpha-beta)*(1:1:T)';
for i = 2:1:T
    sig2t(i) = omega + beta*sig2t(i-1) + alpha*(epst(i-1))^2;
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) APD 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;

