function U_simul = calcU(C_simul,L1_simul,L2_simul,T1_simul,T2_simul,E2_simul,kt_simul,parameters)

eta_cp = parameters.eta_cp ;
MU_0  = parameters.mu_0  ;
psi_l1 = parameters.psi_l1 ;
psi_l2 = parameters.psi_l2 ;
psi_t1 = parameters.psi_t1 ;
psi_t2 = parameters.psi_t2 ;
rhou_l = parameters.rhou_l ;
rhou_t = parameters.rhou_t ;
MU_l1 = parameters.mu_l1 ;
MU_l2 = parameters.mu_l2 ;
MU_t1 = parameters.mu_t1 ;
MU_t2 = parameters.mu_t2 ;
TAU_GAMMA= parameters.tau_gamma;
TT     = parameters.TT      ;
KT     = parameters.KTlarge(1,:);

U_simul=zeros(size(C_simul));
for t=1:TT
	kt   = KT(t);
    mu_0 = MU_0(kt);
    mu_l1= MU_l1(kt);
    mu_l2= MU_l2(kt);
    mu_t1= MU_t1(kt);
    mu_t2= MU_t2(kt);
    tau_gamma = TAU_GAMMA(kt);
    CC_simul = C_simul(:,t) - E2_simul(:,t)*tau_gamma;
    if kt==1
        U_simul(:,t)=mu_0*(1/(1-1/eta_cp))*CC_simul.^(1-1/eta_cp) ... 
        - (1/(1-rhou_l))*(mu_l1*L1_simul(:,t).^(1-1/psi_l1) + mu_l2*L2_simul(:,t).^(1-1/psi_l2)).^(1-rhou_l)...
        - (1/(1-rhou_t))*(mu_t1*T1_simul(:,t).^(1-1/psi_t1) + mu_t2*T2_simul(:,t).^(1-1/psi_t2)).^(1-rhou_t);
    else
        U_simul(:,t)=mu_0*(1/(1-1/eta_cp))*CC_simul.^(1-1/eta_cp) ... 
        - (1/(1-rhou_l))*(mu_l1*L1_simul(:,t).^(1-1/psi_l1) + mu_l2*L2_simul(:,t).^(1-1/psi_l2)).^(1-rhou_l);
    end
end             

end