function  f = SkewTGARCHLogLik(x);
global GARCH_SERIES;
% SKEW_T_GARCH_LOG_LIK.M        Constructs the objective function involved 
%   in the ML estimation of the parameters of a Skew-t GARCH(1,1) 
%   process:
%       Y[t] = mu + sig[t]*eta[t] where
%       sig[t]^2 = alpha0 + alpha1*sig[t-1]^2 + alpha2*(Y[t-1]-mu)^2GARCH(1,1) 
%   with innovations eps[t] that follow a SKEW-T distribution with pdf 
%   (Giot and Laurent 2004):
%           p(x,xi,nu) =
%               (2/(xi+1/xi))*s*c*[1+((xi*(s*x+m))^2)/(nu-2)]^(-1*(nu+1)/2), if x<=-m/s
%           p(x,xi,nu) = 
%               (2/(xi+1/xi))*s*c*[1+(((s*x+m)/xi)^2)/(nu-2)]^(-1*(nu+1)/2), if x>=-m/s
%           with    c = (gamma((nu+1)/2))/(((pi*(nu-2))^0.5)*gamma(nu/2));
%                   m = gamma((nu-1)/2)*((nu-2)^0.5)*(xi-1/xi)/((pi^0.5)*gamma(nu/2));
%                   s = (xi^2+xi^(-2)-1-m^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
xi = exp(x(5,1)); nu = x(6,1);  % skew-t 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) Skew-T log-likelihood for FMINCON
c = (gamma((nu+1)/2))/(((pi*(nu-2))^0.5)*gamma(nu/2));
m = gamma((nu-1)/2)*((nu-2)^0.5)*(xi-1/xi)/((pi^0.5)*gamma(nu/2));
s = (xi^2+xi^(-2)-1-m^2)^0.5;
S1 = ((2/(xi+1/xi))*s*c*[1+((xi*(s*Innovations+m)).^2)/(nu-2)].^(-1*(nu+1)/2)).*(Innovations <= -m/s);
S2 = ((2/(xi+1/xi))*s*c*[1+(((s*Innovations+m)/xi).^2)/(nu-2)].^(-1*(nu+1)/2)).*(Innovations >= -m/s);
S = -1*log(S1+S2);
f = 0.5*sum(log(sig2t))/Tinn + sum(S)/Tinn;

