function [rr, alpha, beta1, beta2,  beta3] = moment_gen_function_gpu(rho1, v1, k1, s1, rho2, v2, k2, s2, c0, c1, c2, muy, sy, mv, rhoj, 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, now hardcoded for gpu to be 2F_Gaussian
%       x:      model parameters, now input individually
%       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
%

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); 

t = a;

w1=0; w2=0; w3=real(phi1); w4=imag(phi1); w5=real(phi2); w6=imag(phi2);

h = (b-a)/N;

for i = 1:N
    [k11,k12,k13,k14,k15,k16] = ODEsystem_2F_Gaussian_gpu(h, t , w1,w2,w3,w4,w5,w6,rho1, v1, k1, s1, rho2, v2, k2, s2, c0, c1, c2, muy, sy, mv, rhoj, uu, rf, d);
    [k21,k22,k23,k24,k25,k26] = ODEsystem_2F_Gaussian_gpu(h, t+h/2, w1+0.5*k11, w2+0.5*k12, w3+0.5*k13, w4+0.5*k14, w5+0.5*k15, w6+0.5*k16,rho1, v1, k1, s1, rho2, v2, k2, s2, c0, c1, c2, muy, sy, mv, rhoj, uu, rf, d);
    [k31,k32,k33,k34,k35,k36] = ODEsystem_2F_Gaussian_gpu(h, t+h/2, w1+0.5*k21, w2+0.5*k22, w3+0.5*k23, w4+0.5*k24, w5+0.5*k25, w6+0.5*k26,rho1, v1, k1, s1, rho2, v2, k2, s2, c0, c1, c2, muy, sy, mv, rhoj, uu, rf, d);
    [k41,k42,k43,k44,k45,k46] = ODEsystem_2F_Gaussian_gpu(h, t+h, w1+k31, w2+k32, w3+k33, w4+k34, w5+k35, w6+k36,rho1, v1, k1, s1, rho2, v2, k2, s2, c0, c1, c2, muy, sy, mv, rhoj, uu, rf, d);
    w1 = w1 + (k11 + 2*k21 + 2*k31 + k41)/6;
    w2 = w2 + (k12 + 2*k22 + 2*k32 + k42)/6;
    w3 = w3 + (k13 + 2*k23 + 2*k33 + k43)/6;
    w4 = w4 + (k14 + 2*k24 + 2*k34 + k44)/6;
    w5 = w5 + (k15 + 2*k25 + 2*k35 + k45)/6;
    w6 = w6 + (k16 + 2*k26 + 2*k36 + k46)/6;
    t = a + i*h;
end

alpha = w1 + 1i*w2;
beta1 = w3 + 1i*w4;
beta2 = w5 + 1i*w6;
beta3 = 0.*w1 + 0.*1i*w2;

rr = exp(alpha + beta1.*v0_1 + beta2.*v0_2);  
