function llf_vec = llf_exceedance_tv_stan(Y_data,nobs_vec,xsi_vec,sig_vec,mu_vec,threshold_vec)
    % Calculate Log-likehood function for PGP model using the approximation from STAN
    % Paramters are time varying
    T = size(Y_data,1);
    llf_vec = zeros(T,1);
    
    for t = 1:T
        xsi = xsi_vec(t);
        sig = sig_vec(t);
        mu = mu_vec(t);
        k = nobs_vec(t);
        llf_t = 0;
        tau=exp(-6.0*abs(xsi));
        % k part of likelihood
        x = threshold_vec(t);
        z=(x-mu)/sig;
        if (1+xsi*z<tau)
		      	lht=-log(tau)/xsi-(z-(tau-1)/xsi)/tau;
		    else
			    lht=-log1p(xsi*z)/xsi;
        end
        llf_t = -exp(lht)-log(factorial(k));
        % other part of likelihood
        for i = 1:k
            x = Y_data(t,i);
            z=(x-mu)/sig;
            if (1+xsi*z<tau)
                lht=-log(tau)/xsi-(z-(tau-1)/xsi)/tau;
           else
                lht=-log1p(xsi*z)/xsi;
           end
           llf_t = llf_t + (1+xsi)*lht-log(sig);
        end
        llf_vec(t) = llf_t;
    end
    
end