% This function takes as input the different simulated types
% and generates a simulated dataset that includes households of different types
function simul_output = simul_all_het(parameters,simul_output_large)

mu_l2_shift   = parameters.mu_l2_shift;
P_mu_l2_shift = parameters.P_mu_l2_shift;
KTlarge       = parameters.KTlarge;
P_KTlarge     = parameters.P_KTlarge;

nk = size(KTlarge,1);
ns = size(mu_l2_shift,2);
n  = size(simul_output_large(1,1).F1_simul,1);

% Number of types is (nk-1)*ns. 
F1_simul = [];
F2_simul = [];
u1_simul = [];
u2_simul = [];
v1_simul = [];
v2_simul = [];
w0_1_simul = [];
w0_2_simul = [];
lw_1_simul = [];
lw_2_simul = [];
Y1_simul = [];
Y2_simul = [];
Y_simul  = [];
atY_simul= [];
C_simul  = [];
A_simul  = [];
H1_simul = [];
H2_simul = [];
L1_simul = [];
L2_simul = [];
T1_simul = [];
T2_simul = [];
kt_simul = [];

for k=1:nk
    for u=1:ns
        nn = ceil(n*P_KTlarge(k)*P_mu_l2_shift(u));
        F1_simul = [F1_simul; simul_output_large(k,u).F1_simul(1:nn,:)]; 
        F2_simul = [F2_simul; simul_output_large(k,u).F2_simul(1:nn,:)]; 
        u1_simul = [u1_simul; simul_output_large(k,u).u1_simul(1:nn,:)]; 
        u2_simul = [u2_simul; simul_output_large(k,u).u2_simul(1:nn,:)]; 
        v1_simul = [v1_simul; simul_output_large(k,u).v1_simul(1:nn,:)]; 
        v2_simul = [v2_simul; simul_output_large(k,u).v2_simul(1:nn,:)]; 
        w0_1_simul = [w0_1_simul; simul_output_large(k,u).w0_1_simul(1:nn,:)]; 
        w0_2_simul = [w0_2_simul; simul_output_large(k,u).w0_2_simul(1:nn,:)]; 
        lw_1_simul = [lw_1_simul; simul_output_large(k,u).lw_1_simul(1:nn,:)]; 
        lw_2_simul = [lw_2_simul; simul_output_large(k,u).lw_2_simul(1:nn,:)]; 
        Y1_simul = [Y1_simul; simul_output_large(k,u).Y1_simul(1:nn,:)]; 
        Y2_simul = [Y2_simul; simul_output_large(k,u).Y2_simul(1:nn,:)]; 
        Y_simul  = [Y_simul; simul_output_large(k,u).Y_simul(1:nn,:)]; 
        atY_simul  = [atY_simul; simul_output_large(k,u).atY_simul(1:nn,:)]; 
        C_simul  = [C_simul; simul_output_large(k,u).C_simul(1:nn,:)]; 
        A_simul  = [A_simul; simul_output_large(k,u).A_simul(1:nn,:)]; 
        H1_simul = [H1_simul; simul_output_large(k,u).H1_simul(1:nn,:)]; 
        H2_simul = [H2_simul; simul_output_large(k,u).H2_simul(1:nn,:)]; 
        L1_simul = [L1_simul; simul_output_large(k,u).L1_simul(1:nn,:)]; 
        L2_simul = [L2_simul; simul_output_large(k,u).L2_simul(1:nn,:)]; 
        T1_simul = [T1_simul; simul_output_large(k,u).T1_simul(1:nn,:)]; 
        T2_simul = [T2_simul; simul_output_large(k,u).T2_simul(1:nn,:)]; 
        kt_simul = [kt_simul; simul_output_large(k,u).kt_simul(1:nn,:)]; 
    end
end

% Case without kids (always last in KT)
F1_simul_nk = [];
F2_simul_nk = [];
u1_simul_nk = [];
u2_simul_nk = [];
v1_simul_nk = [];
v2_simul_nk = [];
w0_1_simul_nk = [];
w0_2_simul_nk = [];
lw_1_simul_nk = [];
lw_2_simul_nk = [];
Y1_simul_nk = [];
Y2_simul_nk = [];
Y_simul_nk  = [];
atY_simul_nk= [];
C_simul_nk  = [];
A_simul_nk  = [];
H1_simul_nk = [];
H2_simul_nk = [];
L1_simul_nk = [];
L2_simul_nk = [];
T1_simul_nk = [];
T2_simul_nk = [];
kt_simul_nk = [];

for u=1:ns
    nn = ceil(n*P_mu_l2_shift(u));
    F1_simul_nk = [F1_simul_nk; simul_output_large(nk,u).F1_simul(1:nn,:)]; 
    F2_simul_nk = [F2_simul_nk; simul_output_large(nk,u).F2_simul(1:nn,:)]; 
    u1_simul_nk = [u1_simul_nk; simul_output_large(nk,u).u1_simul(1:nn,:)]; 
    u2_simul_nk = [u2_simul_nk; simul_output_large(nk,u).u2_simul(1:nn,:)]; 
    v1_simul_nk = [v1_simul_nk; simul_output_large(nk,u).v1_simul(1:nn,:)]; 
    v2_simul_nk = [v2_simul_nk; simul_output_large(nk,u).v2_simul(1:nn,:)]; 
    w0_1_simul_nk = [w0_1_simul_nk; simul_output_large(k,u).w0_1_simul(1:nn,:)]; 
    w0_2_simul_nk = [w0_2_simul_nk; simul_output_large(k,u).w0_2_simul(1:nn,:)]; 
    lw_1_simul_nk = [lw_1_simul_nk; simul_output_large(nk,u).lw_1_simul(1:nn,:)]; 
    lw_2_simul_nk = [lw_2_simul_nk; simul_output_large(nk,u).lw_2_simul(1:nn,:)]; 
    Y1_simul_nk = [Y1_simul_nk; simul_output_large(nk,u).Y1_simul(1:nn,:)]; 
    Y2_simul_nk = [Y2_simul_nk; simul_output_large(nk,u).Y2_simul(1:nn,:)]; 
    Y_simul_nk  = [Y_simul_nk; simul_output_large(nk,u).Y_simul(1:nn,:)]; 
    atY_simul_nk  = [atY_simul_nk; simul_output_large(nk,u).atY_simul(1:nn,:)]; 
    C_simul_nk  = [C_simul_nk; simul_output_large(nk,u).C_simul(1:nn,:)]; 
    A_simul_nk  = [A_simul_nk; simul_output_large(nk,u).A_simul(1:nn,:)]; 
    H1_simul_nk = [H1_simul_nk; simul_output_large(nk,u).H1_simul(1:nn,:)]; 
    H2_simul_nk = [H2_simul_nk; simul_output_large(nk,u).H2_simul(1:nn,:)]; 
    L1_simul_nk = [L1_simul_nk; simul_output_large(nk,u).L1_simul(1:nn,:)]; 
    L2_simul_nk = [L2_simul_nk; simul_output_large(nk,u).L2_simul(1:nn,:)]; 
    T1_simul_nk = [T1_simul_nk; simul_output_large(nk,u).T1_simul(1:nn,:)]; 
    T2_simul_nk = [T2_simul_nk; simul_output_large(nk,u).T2_simul(1:nn,:)]; 
    kt_simul_nk = [kt_simul_nk; simul_output_large(nk,u).kt_simul(1:nn,:)]; 
end


%% Assign to struct
simul_output.F1_simul = F1_simul;
simul_output.F2_simul = F2_simul;
simul_output.u1_simul = u1_simul;
simul_output.u2_simul = u2_simul;
simul_output.v1_simul = v1_simul;
simul_output.v2_simul = v2_simul;
simul_output.w0_1_simul = w0_1_simul;
simul_output.w0_2_simul = w0_2_simul;
simul_output.lw_1_simul = lw_1_simul;
simul_output.lw_2_simul = lw_2_simul;
simul_output.Y1_simul  = Y1_simul;
simul_output.Y2_simul  = Y2_simul;
simul_output.Y_simul  = Y_simul;
simul_output.atY_simul= atY_simul;
simul_output.C_simul  = C_simul;
simul_output.A_simul  = A_simul;
simul_output.H1_simul = H1_simul;
simul_output.H2_simul = H2_simul;
simul_output.L1_simul = L1_simul;
simul_output.L2_simul = L2_simul;
simul_output.T1_simul = T1_simul;
simul_output.T2_simul = T2_simul;
simul_output.kt_simul = kt_simul;

simul_output.F1_simul_nk = F1_simul_nk;
simul_output.F2_simul_nk = F2_simul_nk;
simul_output.u1_simul_nk = u1_simul_nk;
simul_output.u2_simul_nk = u2_simul_nk;
simul_output.v1_simul_nk = v1_simul_nk;
simul_output.v2_simul_nk = v2_simul_nk;
simul_output.w0_1_simul_nk = w0_1_simul_nk;
simul_output.w0_2_simul_nk = w0_2_simul_nk;
simul_output.lw_1_simul_nk = lw_1_simul_nk;
simul_output.lw_2_simul_nk = lw_2_simul_nk;
simul_output.Y1_simul_nk  = Y1_simul_nk;
simul_output.Y2_simul_nk  = Y2_simul_nk;
simul_output.Y_simul_nk  = Y_simul_nk;
simul_output.atY_simul_nk= atY_simul_nk;
simul_output.C_simul_nk  = C_simul_nk;
simul_output.A_simul_nk  = A_simul_nk;
simul_output.H1_simul_nk = H1_simul_nk;
simul_output.H2_simul_nk = H2_simul_nk;
simul_output.L1_simul_nk = L1_simul_nk;
simul_output.L2_simul_nk = L2_simul_nk;
simul_output.T1_simul_nk = T1_simul_nk;
simul_output.T2_simul_nk = T2_simul_nk;
simul_output.kt_simul_nk = kt_simul_nk;

