function phi = CharFun(params, r, tau, u, model, JT)
% Characteristic function of the jump component 
if JT == 1      % Merton
    % Parameteres: lambda, muj, sigmaj
    lambda = params(6);
    muj = params(7);
    sigmaj = params(8);    
    psi_d = @(u) lambda*(exp(1i*u*muj)-0.5*u.^2*sigmaj^2-1);
else
    % Parameteres: lambda, psip, psim
    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
        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);
end

muj_t = -psi_d(-1i);    % risk-neutral compensator
phi_d = exp(tau*(psi_d(u)+1i.*u.*muj_t));

% Characteristic function of the diffusive component
% Parameters: rho, beta, sigma, alpha, eta
rho = params(1);
beta = params(2);
sigma = params(3);
alpha = params(4);
eta = params(5);

v = u.*sigma*sqrt(tau);
mu_t = r-sigma^2/2;

if strcmp(model,'Edgeworth')
    % Local Edgeworth expansion
    first_adj = -1i*v.^3*beta*rho*sqrt(tau)/(2*sigma);
    second_adj = -v.^2*(alpha/(2*sigma)+beta.^2/(4*sigma^2))*tau;
    third_adj = beta^2/(24*sigma^2)*v.^2.*(4*v.^2-rho^2*v.^2.*(3*v.^2-8))*tau;
    fourth_adj = eta/(6*sigma)*v.^4*tau; 
    phi_c = exp(1i*v*mu_t*tau./(sigma*sqrt(tau))-v.^2/2).*(1+first_adj+second_adj+third_adj+fourth_adj);
elseif strcmp(model,'Gram-Charlier')
    % Local Gram-Charlier expansion
    a = -alpha/(2*sigma)-beta^2/(4*sigma^2);
    b = -beta*rho/(2*sigma);
    c = beta^2/(6*sigma^2)+beta^2*rho^2/(3*sigma^2)+eta/(6*sigma);
    d = -beta^2*rho^2/(8*sigma^2);

    % centered moments
    second_moment = 1-2*a*tau;
    third_moment = -6*b*sqrt(tau);
    fourth_moment = 3*(1-4*a*tau+8*c*tau);
    fifth_moment = -60*b*sqrt(tau);
    sixth_moment = 15*(1-6*a*tau+24*c*tau-48*d*tau);

    % coefficients c_j
    c2 = second_moment-1;
    c3 = third_moment;
    c4 = fourth_moment-6*second_moment+3;
    c5 = fifth_moment-10*third_moment;
    c6 = sixth_moment-15*fourth_moment+45*second_moment-15;
    phi_c = exp(1i*v*mu_t*tau./(sigma*sqrt(tau))-v.^2/2).*(1+c2/factorial(2)*(1i.*v).^2+c3/factorial(3)*(1i.*v).^3 + ...
         c4/factorial(4)*(1i.*v).^4 + c5/factorial(5)*(1i.*v).^5 + c6/factorial(6)*(1i.*v).^6);
elseif strcmp(model,'BSM') 
    % Gaussian characteristic function
    phi_c = exp(1i*v*mu_t*tau./(sigma*sqrt(tau))-v.^2/2);
end

% Complete characteristic function
phi = phi_d.*phi_c;

end