function [llf, h] = loglike_skt_FS(est, r)
% [llf, h] = loglike_skt_FS(est, r)  
% Calculates the loglikelihood assuming it follows a skew t distribution Fernandez Steel (1998)
%
% Author: Christian Contino
% Date:   12 August 2014

    alpha0        = est(1);               % Sets alpha0 to first column value from estimation matrix
    alpha1        = est(2);               % Sets alpha1 to second column value from estimation matrix
    beta1         = est(3);               % Sets beta1 to third column value from estimation matrix
    nu            = est(4);               % Sets nu to fourth column value from estimation matrix
    lambda        = est(5);               % Sets lambda to fifth column from estimation matrix
    mean_ret      = est(6);               % Sets mean of returns
    start_vol     = est(7);               % Day 1 volatility for likelihood
    
    N             = length(r);            % Set loop length N to size of return series a
    h             = zeros(N,1);           % Fills volatility estimates vector with zeros
    a             = r - mean_ret;
    a2            = a.^2;
   
   % Parameter Restrictions         
    if alpha0 < 0
        llf = -inf;
        return;
    end
    
    if alpha1 <= 0
        llf = -inf;
        return;
    end
    
    if beta1 <= 0
        llf = -inf;
        return;
    end
    
    if alpha1 + beta1 >= 1
        llf = -inf;
        return;
    end
    
    if nu <= 2
        llf = -inf;
        return;
    end
    
    if nu > 200
        llf = -inf;
        return;
    end
        
    if lambda < 0
        llf = -inf;
        return;
    end
    
    if start_vol >= 200
        llf = -inf;
        return
    end
    
    if start_vol <= 0
        llf = -inf;
        return
    end   
   
    for t = 1:N
        if (t == 1)                                     % If T is equal to 1. Needed as we don't have previous observations
            h(t) = start_vol;                        % Sets volatility ln_h to start_vol parameter
        else
            h(t) = alpha0 + alpha1 .* a2(t-1) + beta1 .* h(t-1); % Calculates volatility h based on proposed parameters and squared returns a2 using GARCH(1,1) equation
        end
    end 

m =  ((gamma((nu-1)/2)*sqrt(nu-2))./(gamma(nu/2)*sqrt(pi))) .* (lambda - 1./lambda);
l2  = lambda.^2;
s   = sqrt(l2 + 1./l2 - 1 - m.^2);

find1 = (a>=(-m./s));
find2 = (a<(-m./s));

prior = log(1/nu); %Inverse Prior on DoF
LL1 = N.*(log(gamma((nu+1)./2)) - log(gamma(nu./2)) - 0.5.*log((nu-2).*pi) + log(2.*lambda.*s) - log(1+l2));
LL2 = (((nu+1)./2) .* log(1 + (lambda.^(-2).*(m.*sqrt(h)+s.*a).^2)./ ((nu-2).*h))); 
LL3 = (((nu+1)./2) .* log(1 + (lambda.^(2) .*(m.*sqrt(h)+s.*a).^2)./ ((nu-2).*h))); 
llf = LL1 - sum(LL2.*(find1))  - sum(LL3.*(find2)) - 0.5.*sum(log(h)) + prior;      
end