function [c1, c2, c4] = GetCumulants(r, params, T, JT)
% Computation of the cumulants for different jump types
if JT == 1  % Merton
    sigma = params(3);
    lambda = params(6);
    muj = params(7);
    sigmaj = params(8);
    psi_d = @(u) lambda*(exp(1i*u*muj)-0.5*u.^2*sigmaj^2-1);
    rn_drift = r-sigma^2/2-psi_d(-1i);

    c1 = (rn_drift+lambda*muj)*T;
    c2 = (sigma^2+lambda*sigmaj^2+lambda*muj^2)*T;
    c4 = lambda*(3*sigmaj^3+6*muj^2*sigmaj^2+muj^4)*T;
else    
    sigma = params(3);
    lambda = params(6);
    psip = params(7);
    psim = params(8);
    if JT == 2      % Kou
        nu = -1;
        psi_dm = @(u) gamma(-nu).*psim.^nu.*lambda.*((1+1i.*u/psim).^nu-1-1i.*u*nu/psim);
        psi_dp = @(u) gamma(-nu).*psip.^nu.*lambda.*((1-1i.*u/psip).^nu-1+1i.*u*nu/psip);
    elseif JT == 3      % VG
        nu = 0;
        psi_dm = @(u) -lambda*(-1i.*u/psim+log(1+1i.*u/psim));
        psi_dp = @(u) -lambda*(1i.*u/psip+log(1-1i.*u/psip));
    else    % NIG
        nu = 1/2;
        psi_dm = @(u) gamma(-nu).*psim.^nu.*lambda.*((1+1i.*u/psim).^nu-1-1i.*u*nu/psim);
        psi_dp = @(u) gamma(-nu).*psip.^nu.*lambda.*((1-1i.*u/psip).^nu-1+1i.*u*nu/psip);
    end
    psi_d = @(u) psi_dm(u)+psi_dp(u);
    rn_drift = r-sigma^2/2-psi_d(-1i);
    c1 = rn_drift*T;
    c2 = T*(sigma^2+gamma(2-nu)*lambda*psip^(nu-2)+gamma(2-nu)*lambda*psim^(nu-2));
    c4 = T*(gamma(4-nu)*lambda*psip^(nu-4)+gamma(4-nu)*lambda*psim^(nu-4));
end

end