function [r_real, r_imag] = DE_MGF(u1, u2, u3, flag_mp, lambdam, lambdap, muv1_m, muv1_p, muv3_m, muv3_p, ls_flag)
% Moment generating function for double-exponential jump component
%
% Computes real and imaginary parts of the MGF contribution from positive
% or negative jumps with optional leverage (ls_flag), uses the Faddeeva
% function W.m for the leverage case

r_real = NaN;
r_imag = NaN;

%% Negative Jumps

if flag_mp == -1 
    if  lambdam == 0 && muv1_m == 0
        r = 1;
    elseif (u2 == 0 || muv1_m == 0) && (u3 == 0 || muv3_m == 0)
        r = lambdam./(lambdam+u1);
    else
        if ls_flag == 1
            am = -(u2 * muv1_m + u3 * muv3_m);
            bm = (u1 + lambdam)./(2*(u2 * muv1_m + u3 * muv3_m)); 
            % Exp[z^2] (Int_{-Inf}^{z} Exp[-y^2] dy) =
            % Exp[z^2] (sqrt(pi)/2 + erf(z)*sqrt(pi)/2) = 
            % Exp[z^2]  sqrt(pi)/2 (1+erf(z)) =
            % Exp[z^2]  sqrt(pi)/2 (1-erf(-z)) =
            % Exp[z^2]  sqrt(pi)/2 erfc(-z) =
            % sqrt(pi)/2 W(-iz)
            r = lambdam*(1/sqrt(am)) * W(-1i*sqrt(am)*bm) * sqrt(pi)/2;
        elseif ls_flag  == 0
            r = lambdam./(u1 - u2.*muv1_m + lambdam);
        end
    end
end

%% Positive Jumps

if flag_mp == 1
    if lambdap == 0 && muv1_p == 0
        r = 1;
    elseif (u2 == 0 || muv1_p == 0) && (u3 == 0 || muv3_p == 0)
        r = lambdap./(lambdap-u1);
    else
        if ls_flag == 1
            bp = (u1 - lambdap)./(2*(u2 * muv1_p + u3 * muv3_p)); 
            ap = -(u2 * muv1_p + u3 * muv3_p);
            % Exp[z^2] (Int_{z}^{+Inf} Exp[-y^2] dy) =
            % Exp[z^2]  erfc(z)*sqrt(pi)/2) = 
            % sqrt(pi)/2 W(iz)
            r = lambdap*(1/sqrt(ap)) * W(1i*sqrt(ap)*bp) * sqrt(pi)/2;
        elseif ls_flag  == 0
            r = -lambdap./(u1 + u2.*muv1_p - lambdap);
        end
    end
end

r_real = real(r);
r_imag = imag(r);

