function [rr, alpha, beta1, beta2,  beta3] = moment_gen_function(distr, x, uu, phi1, phi2, phi3, sv1, sv2, sv3, T, rf, d)
%
% moment_gen_function.m computes the Joint Conditional MGF (or the
% Characteristic Function -CF-) of log(S_T/S_t), V1_T, V2_T, and V3_T
% corresponding to the model "distr" at points uu, phi1, phi2, and phi3. 
%
% Note:
% CF_t(x,u,phi1,phi2,phi3,T) =
% E_t[exp(uu*y_T+phi1*V_{1,T}+phi2*V_{2,T}+phi3*V_{3,T})] =
% exp(alpha(x,T) + beta1(x,T)V_{1,t} + beta1(x,T)V_{2,t} + beta3(x,T)V_{3,t})
% OUTPUT: 
%
%           M(uu,phi1,phi2,t,T) = E_t[exp(uu*y_T+ alpha + phi1*V_{1,T} + phi2*V_{2,T})] 
% INPUT:
%       distr:  model distribution
%       x:      model parameters
%       uu:     evalutation points for y_T (can be a vector)
%       phi1:   evalutation points for V_1 (can be a vector)
%       phi2:   evalutation points for V_2 (can be a vector)
%       phi3:   evalutation points for V_3 (can be a vector)
%       sv1:    initial value of the first factor
%       sv2:    initial value of the second factor
%       sv3:    initial value of the third factor
%       T:      maturity (in years, convention 365 days per year)
%       rf:     interest rate
%       d:      dividend yield
%
% Nicola Fusari (2015):
% You can distribute and use the code. Please cite the author if you use
% the code in in your research.
% 
% The code used to Solve the Systme of ODEs is from Alain Kapitho
% (University of Pretoria), available at:
% http://m2matlabdb.ma.tum.de/rk4_systems.m?MP_ID=413
%

rr = zeros(1, length(uu));

ndays = round(T*365);

TT = (0:ndays)'/365;

alpha = complex(zeros(ndays+1,length(uu)));
beta1 = complex(zeros(ndays+1,length(uu)));
beta2 = complex(zeros(ndays+1,length(uu)));
beta3 = complex(zeros(ndays+1,length(uu)));

switch distr 
    case 'BS'
      
        alpha = uu.*(rf-d).*T;
        beta1 = (-(uu/2)+(uu.^2)./2).*T; 
        
        rr = exp(alpha(end,:)+beta1(end,:).*sv1);
        
    case 'Heston_Merton'
        
        % In fraction of days
        dt = T;

        rho       = x(1);
        v         = x(2);
        kappa_v   = x(3);
        sigma_v   = x(4);
        lambda_y0 = x(9);
        lambda_y1 = x(10);
        mu_y      = x(12);   
        sigma_y   = x(13);  

        % Compute the jumps compensators:
        m_y = exp(mu_y+sigma_y^2*0.5) - 1;

        if (kappa_v == 0)  && (lambda_y1==0)
             alpha = uu.*((rf-d)- 0.5.*v).*dt + 0.5.*uu.^2.*v.*dt + lambda_y0.*dt.*(exp(mu_y.*uu+sigma_y^2*0.5.*uu.^2) - 1 -uu.*m_y);
            rr = exp(alpha(end,:));
        else

            b = sigma_v.*rho.*uu - kappa_v;
            a = uu.*(1-uu) -2.*lambda_y1.*(exp(uu.*mu_y+0.5.*uu.^2.*sigma_y.^2)-1-uu.*m_y);

            zz = b.^2+a.*sigma_v^2;
            g = sqrt(abs(zz)) .*exp(1i.*atan2(imag(zz),real(zz))./2);

            beta1 = -(a.*(1-exp(-g.*dt)))./(2*g-(g+b).*(1-exp(-g.*dt)));

            z = 1-((g+b)./(2*g)).*(1-exp(-g.*dt)); 
            logz = log(abs(z)) + 1i.* atan2(imag(z),real(z));

            alpha_0 = (rf-d).*uu.*dt - v.*kappa_v.*(  dt.*(g+b)./(sigma_v.^2) + 2.*logz./sigma_v.^2 );

            alpha = alpha_0 + lambda_y0.*dt.*(exp(mu_y.*uu+(sigma_y.^2.*uu.^2)./2)-1-m_y.*uu);

            rr = exp(alpha(end,:)+beta1(end,:).*sv1);

        end
        
    case 'DPS'
     
        
        % In fraction of days
        dt = T;

        rho       = x(1);
        v         = x(2);
        kappa_v   = x(3);
        sigma_v   = x(4);
        lambda_y0 = x(9);
        mu_y      = x(12);   
        sigma_y   = x(13);  
        mu_v      = x(14);
        rho_j     = x(15);

        
        % Compute the correlated jumps compensator:
        comp1 = (exp(mu_y+0.5*sigma_y^2)/(1-rho_j*mu_v)) - 1.0;

        b = sigma_v.*rho.*uu - kappa_v;
        a = uu.*(1-uu);
        c = 1-rho_j.*mu_v.*uu;

        zz = b.^2+a.*sigma_v^2;
        g = sqrt(abs(zz)) .*exp(1i.*atan2(imag(zz),real(zz))./2);
        
        
        LD = 1 - (((g+b).*c-mu_v.*a)./(2.*g.*c)).*(1-exp(-g.*dt));
        logD = log(abs(LD)) + 1i.* atan2(imag(LD),real(LD));

        
        D = complex(zeros(1,length(uu)));

        if lambda_y0>0 && mu_v>0
            for h=1:length(uu)
                if uu(h)==0
                    D(h) = dt.*((g(h)-b(h))./((g(h)-b(h)).*c(h)+mu_v.*a(h)));
                else
                    D(h) = dt.*((g(h)-b(h))./((g(h)-b(h)).*c(h)+mu_v.*a(h)))  - ((2.*mu_v.*a(h))./((g(h).*c(h)).^2-(b(h).*c(h)-mu_v.*a(h)).^2)).*logD(h) ;
                end
            end
        elseif lambda_y0>0 && mu_v==0
            D = dt ;
        elseif lambda_y0==0
            D = 0;
        end

        beta1 = -(a.*(1-exp(-g.*dt)))./(2*g-(g+b).*(1-exp(-g.*dt)));

        z = 1-((g+b)./(2*g)).*(1-exp(-g.*dt)); 
        logz = log(abs(z)) + 1i.* atan2(imag(z),real(z));

        alpha_0 = -rf*dt + (rf-d).*uu.*dt - v.*kappa_v.*(  dt.*(g+b)./(sigma_v.^2) + 2.*logz./(sigma_v.^2) );

        alpha = alpha_0 - lambda_y0.*dt.*(1 +comp1.*uu) +  lambda_y0.*( (exp(mu_y.*uu+(sigma_y.^2.*uu.^2)./2).*D));

        rr = exp(alpha(end,:)+beta1(end,:).*sv1);
    
    case 'Heston_GTS'
        
        % In fraction of days
        dt = T;
        
        % Volatility Factor
        rho     = x(1);
        v       = x(2);
        kappa_v = x(3);
        sigma_v = x(4);

        % Jump Specification
        c0_p = x(5); 
        c0_m = x(6);
        c1_p = x(7);
        c1_m = x(8);
        alpha_p  = x(9);
        alpha_m  = x(10);
        lambda_p = x(11);
        lambda_m = x(12);

        if abs(alpha_p-0)<eps 
            C0_p =    -c0_p.*(uu./lambda_p +log(1-uu./lambda_p));
            Corr0_p = -c0_p.*( 1./lambda_p +log(1- 1./lambda_p));
        elseif abs(alpha_p - 1)<eps
            C0_p =    c0_p.* uu + c0_p.*(lambda_p- uu).*log(1  -  uu/lambda_p);
            Corr0_p = c0_p.* 1.0+ c0_p.*(lambda_p-1.0).*log(1.0- 1.0/lambda_p);
        else
            C0_p =    gamma(-alpha_p).*(lambda_p.^alpha_p) .* c0_p.*( power(complex(real(1-(uu./lambda_p)),imag(1-(uu./lambda_p))),alpha_p) - 1 + (uu.*alpha_p)./(lambda_p));  
            Corr0_p = gamma(-alpha_p).*(lambda_p.^alpha_p) .* c0_p.*( power(complex(real(1-( 1./lambda_p)),imag(1-( 1./lambda_p))),alpha_p) - 1 + ( 1.*alpha_p)./(lambda_p));
        end   


        if abs(alpha_m-0)<eps  
            C0_m =    -c0_m.*(-uu./lambda_m + log(1+uu./lambda_m));
            Corr0_m = -c0_m.*(- 1./lambda_m + log(1+ 1./lambda_m));
        elseif abs(alpha_m - 1)<eps
            C0_m =    -c0_m.*uu + c0_m.*(lambda_m+ uu).*log(1.0+ uu/lambda_m);
            Corr0_m = -c0_m.*1.0+ c0_m.*(lambda_m+1.0).*log(1.0+1.0/lambda_m);
        else
             C0_m =    gamma(-alpha_m).*(lambda_m.^alpha_m) .* c0_m.*( power(complex(real(1+uu./lambda_m),imag(1+uu./lambda_m)),alpha_m) - 1 - (uu.*alpha_m)./(lambda_m)); 
            Corr0_m = gamma(-alpha_m).*(lambda_m.^alpha_m) .* c0_m.*( power(complex(real(1+ 1./lambda_m),imag(1+ 1./lambda_m)),alpha_m) - 1 - ( 1.*alpha_m)./(lambda_m));  
        end

        Corr0 = Corr0_p + Corr0_m;
        C0 = C0_p + C0_m;

        if abs(kappa_v-0)<eps 

            alpha = uu.*(rf-d-Corr0-0.5.*v).*dt + 0.5.*uu.^2.*v.*dt + dt.*C0;
            
            rr = exp(alpha(end,:));
        else

            if abs(alpha_p-0)<eps 
                C1_p =    -c1_p.*(uu./lambda_p +log(1-uu./lambda_p));
                Corr1_p = -c1_p.*( 1./lambda_p +log(1- 1./lambda_p));
            elseif abs(alpha_m -1)<eps
                C1_p =    c1_p.* uu+ c1_p.*(lambda_p- uu).*log(1.0- uu/lambda_p);
                Corr1_p = c1_p.*1.0+ c1_p.*(lambda_p-1.0).*log(1.0-1.0/lambda_p);
            else
                C1_p =    gamma(-alpha_p).*lambda_p.^alpha_p .* c1_p.*((1-uu/lambda_p).^alpha_p - 1+(uu.*alpha_p)./(lambda_p));
                Corr1_p = gamma(-alpha_p).*lambda_p.^alpha_p .* c1_p.*((1- 1/lambda_p).^alpha_p - 1+( 1.*alpha_p)./(lambda_p));
            end   


            if abs(alpha_m-0)<eps  
                C1_m =    -c1_m.*(-uu./lambda_m + log(1+uu./lambda_m));
                Corr1_m = -c1_m.*(- 1./lambda_m + log(1+ 1./lambda_m));
            elseif abs(alpha_m-1)<eps 
                C1_m =    -c1_m.* uu+ c1_m.*(lambda_m+ uu).*log(1.0+ uu./lambda_m);
                Corr1_m = -c1_m.*1.0+ c1_m.*(lambda_m+1.0).*log(1.0+1.0/lambda_m);
            else
                C1_m =    gamma(-alpha_m).*lambda_m.^alpha_m .* c1_m.*((1+uu./lambda_m).^alpha_m - 1 - (uu.*alpha_m)./(lambda_m));  
                Corr1_m = gamma(-alpha_m).*lambda_m.^alpha_m .* c1_m.*((1+ 1./lambda_m).^alpha_m - 1 - ( 1.*alpha_m)./(lambda_m));  
            end

            C1 = C1_p + C1_m;
            Corr1 = Corr1_p + Corr1_m;

            a = uu.*(1-uu) -2.*(C1-Corr1.*uu);

            b = sigma_v.*rho.*uu - kappa_v;
            
            zz = b.^2+a.*sigma_v^2;
            g = sqrt(abs(zz)) .*exp(1i.*atan2(imag(zz),real(zz))./2);

            beta1 = -(a.*(1-exp(-g.*dt)))./(2*g-(g+b).*(1-exp(-g.*dt)));

            alpha = uu.*(rf-d-(Corr0) ).*dt + dt.*C0 - kappa_v.*v.*(dt.*(b+g)./(sigma_v.^2)  +2.*log(1-(g+b)./(2.*g) .*(1-exp(-g.*dt)))./sigma_v^2 );

            rr = exp(alpha(end,:) + beta1(end,:).*sv1);
           
        end
  
    case '2F_Gaussian'
        
        v0_1 = sv1;   
        v0_2 = sv2;   
       
        % The system of ODEs is solved between [a, b]
        a = 0;
        b = T;
        
        % Number of days -> Also, number of steps used to solve the ODEs system:
        N = round(T*365); 

        % Time steps
        t = zeros(1,N+1);
        
        % Matrix that stores the entire solution of the ODEs system for a
        % particular value of uu: 
        w = zeros(6,N+1);

        if N==0
            for g =1:length(uu)
                 beta1(:,g) = real(phi1(g)) + 1i*imag(phi1(g)) ;
                 beta2(:,g) = real(phi2(g)) + 1i*imag(phi2(g));
                
            end
        else
            for g =1:length(uu)

                % Vector of Initial Conditions
                R0 = [0.0; 0.0;  real(phi1(g)) ; imag(phi1(g)) ; real(phi2(g)) ; imag(phi2(g))];

                h = (b-a)/N;        %the step size
                t(1) = a;
                w(:,1) = R0;        %initial conditions

                for i = 1:N
                   k1 = h*ODEsystem_2F_Gaussian(t(i)    , w(:,i)        , x, uu(g), rf, d);
                   k2 = h*ODEsystem_2F_Gaussian(t(i)+h/2, w(:,i)+0.5*k1 , x, uu(g), rf, d);  
                   k3 = h*ODEsystem_2F_Gaussian(t(i)+h/2, w(:,i)+0.5*k2 , x, uu(g), rf, d);  
                   k4 = h*ODEsystem_2F_Gaussian(t(i)+h  , w(:,i)+k3     , x, uu(g), rf, d);  
                   w(:,i+1) = w(:,i) + (k1 + 2*k2 + 2*k3 + k4)/6;
                   t(i+1) = a + i*h;
                end

                % Time x uu
                % each column of alpha, beta1, beta2, and beta3 contain 
                % the solution of the ODE system for a specific value of uu
                alpha(:,g) = w(1,:)' + 1i*w(2,:)';
                beta1(:,g) = w(3,:)' + 1i*w(4,:)';
                beta2(:,g) = w(5,:)' + 1i*w(6,:)';
                beta3(:,g) = 0.*w(1,:)' + 0.*1i*w(2,:)';

            end
        end

        rr = exp(alpha(end,:) + beta1(end,:).*v0_1 + beta2(end,:).*v0_2);  
        
        
    case '3F_DE'
        
        
        v0_1 = sv1;   
        v0_2 = sv2;   
        v0_3 = sv3;   

        % The system of ODEs is solved between [a, b]
        a = 0;
        b = T;
        
        % Number of days -> Also, number of steps used to solve the ODEs system:
        N = round(T*365); 

        % Time steps
        t = zeros(1,N+1);
        
        % Matrix that stores the entire solution of the ODEs system for a
        % particular value of uu: 
        w = zeros(8,N+1);

        for g =1:length(uu)
            
            % Vector of Initial Conditions
            R0 = [0.0; 0.0;  real(phi1(g)) ; imag(phi1(g)) ; real(phi2(g)) ; imag(phi2(g)) ; real(phi3(g)) ; imag(phi3(g))];

            h = (b-a)/N;        %the step size
            t(1) = a;
            w(:,1) = R0;        %initial conditions

            for i = 1:N
               k1 = h*ODEsystem_3F_DE(t(i)    , w(:,i)        , x, uu(g), rf, d);
               k2 = h*ODEsystem_3F_DE(t(i)+h/2, w(:,i)+0.5*k1 , x, uu(g), rf, d);  
               k3 = h*ODEsystem_3F_DE(t(i)+h/2, w(:,i)+0.5*k2 , x, uu(g), rf, d);  
               k4 = h*ODEsystem_3F_DE(t(i)+h  , w(:,i)+k3     , x, uu(g), rf, d);  
               w(:,i+1) = w(:,i) + (k1 + 2*k2 + 2*k3 + k4)/6;
               t(i+1) = a + i*h;
            end

            % Time x uu
            % each column of alpha, beta1, beta2, and beta3 contain 
            % the solution of the ODE system for a specific value of uu
            if imag(uu(g))==0 && imag(phi1(g))==0 && imag(phi2(g))==0 && imag(phi3(g))==0
                alpha(:,g) = w(1,:)';
                beta1(:,g) = w(3,:)';
                beta2(:,g) = w(5,:)';
                beta3(:,g) = w(7,:)';
            else
                alpha(:,g) = w(1,:)' + 1i*w(2,:)';
                beta1(:,g) = w(3,:)' + 1i*w(4,:)';
                beta2(:,g) = w(5,:)' + 1i*w(6,:)';
                beta3(:,g) = w(7,:)' + 1i*w(8,:)';

            end
        end

        rr = exp(alpha(end,:) + beta1(end,:).*v0_1 + beta2(end,:).*v0_2 + beta3(end,:).*v0_3);  
     
    
    case 'HN'
        if length(x) == 4
            
            % Heston & Nandi (2000):  
            % x(1) w, 
            % x(2) b, 
            % x(3) a
            % x(4) c, 

            v0_1 = sv1/365;  
            
           
            for t = ndays:-1:1
                alpha(t,:) = alpha(t+1,:) + uu.*(rf-d) + beta1(t+1,:).*x(1) - 0.5.*log(1-2.*x(3).*beta1(t+1,:));
                beta1(t,:) = uu.*(-0.5 + x(4)) - 0.5.*(x(4)^2) + ( ((uu-x(4)).^2)./(2*(1-2.*x(3).*beta1(t+1,:))) ) + x(2).*beta1(t+1,:);
            end

            rr = exp(alpha(1,:) + beta1(1,:).*v0_1);  
        else
            fprintf('You did not provide the rigth number of parameters: 8');
        end
        
        
    case 'CGARCH'
        if length(x) == 7
                  
            % Christoffersen et al. (2008):  
            % x(1) bs, 
            % x(2) as, 
            % x(3) cs, 
            % x(4) w, 
            % x(5) bl, 
            % x(6) al, 
            % x(7) cl,
            

            v0_1 = sv1/365;
            v0_2 = sv2/365;
            
            for t = ndays:-1:1
                alpha(t,:) = alpha(t+1,:) + (rf-d).*uu - (x(2).*beta1(t+1,:) + x(6).*beta2(t+1,:)) -0.5.*log(1-2.*x(2).*beta1(t+1,:)-2.*x(6).*beta2(t+1,:)) + beta2(t+1,:).*x(4);
                beta1(t,:) = beta1(t+1,:).* (x(1)+x(2)*x(3)^2) - 0.5.*uu + 2.*( ((x(2).*x(3).*beta1(t+1,:) + x(6).*x(7).*beta2(t+1,:) - 0.5.*uu).^2)./(1-2.*x(2).*beta1(t+1,:)-2.*x(6).*beta2(t+1,:)));
                beta2(t,:) = beta2(t+1,:).*x(5) - 0.5.*uu + 2.*( ((x(2).*x(3).*beta1(t+1,:) + x(6).*x(7).*beta2(t+1,:) - 0.5.*uu).^2)./(1-2.*x(2).*beta1(t+1,:)-2.*x(6).*beta2(t+1,:)));
            end
               
            beta2 = beta2-beta1;
            rr = exp(alpha(1,:) + beta1(1,:).*v0_1 + beta2(1,:).*v0_2);  
        else
            fprintf('You did not provide the rigth number of parameters: 12');
        end
            
    
    case 'LHARG'
        
        alpha = complex(zeros(ndays+1,length(uu)));
        beta1 = complex(zeros(ndays+1,length(uu)));
        beta2 = complex(zeros(ndays+1,length(uu)));
        beta3 = complex(zeros(ndays+1,length(uu)));
       
        if length(x) == 10
                 
            % x(1)  d
            % x(2)  delta (shape)
            % x(3)  theta (scale)
            % x(4)  bd
            % x(5)  bw
            % x(6)  bm
            % x(7)  ad
            % x(8)  aw
            % x(9)  am
            % x(10) gamma (laverage)
            
            d     = x(1); 
            del   = x(2); % (shape)
            thet  = x(3); % (scale)
            bd    = x(4);
            bw    = x(5); 
            bm    = x(6);
            ad    = x(7);
            aw    = x(8);
            am    = x(9);
            gm   = x(10);

            RV = sv1/365;
            Y  = sv2/365;
            
            % Reshape \beta and \alpha parameters
            bb = zeros(1,22);
            bb(1)    = bd;
            bb(2:5)  = bw/4;
            bb(6:22) = bm/17;
            
            aa = zeros(1,22);
            aa(1)    = ad;
            aa(2:5)  = aw/4;
            aa(6:22) = am/17;
            
           
            Ld = max((Y./sqrt(RV) - gm.*sqrt(RV)).^2,0);
            
            A = complex(zeros(ndays+1,length(uu), 1));
            X = complex(zeros(ndays+1,length(uu), 1));
            B = complex(zeros(ndays+1,length(uu),22));
            C = complex(zeros(ndays+1,length(uu),22));
            
            lmbd = -0.5;
            
            
            for t = ndays:-1:1
                X(t+1,:) = (uu.*(lmbd) + B(t+1,:,1) + (0.5.*(uu.^2)+gm^2.*C(t+1,:,1)-2.*C(t+1,:,1).*gm.*uu)./(1-2.*C(t+1,:,1)));
        

                A(t,:) = A(t+1,:) + uu.*rf - 0.5*log(1-2*C(t+1,:,1)) - del*f_W(X(t+1,:),thet);
                 
                for n=1:22
                    if n<22
                        B(t,:,n) = B(t+1,:,n+1)+f_V(X(t+1,:),thet).*bb(n);
                        C(t,:,n) = C(t+1,:,n+1)+f_V(X(t+1,:),thet).*aa(n);
                    elseif n==22
                        B(t,:,n) =              f_V(X(t+1,:),thet).*bb(n);
                        C(t,:,n) =              f_V(X(t+1,:),thet).*aa(n);
                    end
                end


            end

            if length(uu)==1
                rr = exp(A(1,:) + flip(RV)'*squeeze(B(1,:,:))  + flip(Ld)'*squeeze(C(1,:,:))); 
            else
                rr = exp(A(1,:) + transpose(sum(repmat(flip(RV)',length(uu),1).*squeeze(B(1,:,:)),2)) + transpose(sum(repmat(flip(Ld)',length(uu),1).*squeeze(C(1,:,:)),2)));
            end

        else
            fprintf('You did not provide the rigth number of parameters: 12');
        end
            
    
    otherwise    
        fprintf('You did not choose any allowed Model')
        
end



%% Helper functions for HARG model

function r = f_V(x,thet)

r = (thet.*x)./(1-thet.*x);


function r = f_W(x,thet)

r = log(1-x.*thet);



