% The path of residual file must be specified ..

% EMPLOC Function evaluation
function [ZZ]= NLLS_EMPLOC_SIM(y,EMP_lag_1,EMPactual_1,sumHHzone,ZoneSize,c,dimZone,dimHH,dimEMP);

%Initialization
lambda=y(1,1); omega=y(2,1); rho=y(3,1); 
aa=y(4,1); bb=y(5,1);

k=1;
dimEMP=1;

%Calculate A
Ainv=zeros(dimZone,dimEMP);
for i=1:dimZone;
        for j=1:dimZone;
            Ainv(i,k)=Ainv(i,k)+(EMP_lag_1(j,k).^aa(1,k).*ZoneSize(j,1).^bb(1,k).*c(i,j).^omega(1,k).*exp(rho(1,k)*c(i,j)));
        end;
        if(Ainv(i,k) < 0.001) Ainv(i,k) = 0.001; end;
        A(i,k)=1./Ainv(i,k);
end;

%Calculate W_EMPLOC
for j=1:dimZone;
        W_EMPLOC(j,k)= EMP_lag_1(j,k).^aa(1,k).*ZoneSize(j,1).^bb(1,k);     
end;

%Calculate modeled EMP
EMPsum=zeros(dimZone,dimEMP);
EMP = zeros(dimZone,dimEMP);
for j=1:dimZone;
        for i=1:dimZone;
            EMPsum(j,k)=EMPsum(j,k)+sumHHzone(i,1).*A(i,k).*W_EMPLOC(j,k).*c(i,j).^omega(1,k).*exp(rho(1,k).*c(i,j));
        end;
        EMP(j,k)=(lambda(1,k).*EMPsum(j,k))+((1-lambda(1,k))*EMP_lag_1(j,k));
        if(EMP(j,k) < 0.001)
            EMP(j,k) = 0.001;
        end;
end;

% Scaling 
sum_EMPactual = sum(EMPactual_1,1);
sum_EMP = sum(EMP,1);
EMP = EMP.*(sum_EMPactual/sum_EMP);


%Calculate residual
k = 1;
epsilon=zeros(dimZone,1);
for i=1:dimZone;
        epsilon(i,k)=EMPactual_1(i,k) - EMP(i,k);
end;
csvwrite('C:\G-LUM Finale\Calibration\Output\epsilon_emp.dat',epsilon);% Output path must be specified

% Likelihood functions ......
L = 0;
for i=1:dimZone;
    L = L + (EMPactual_1(i,k)*log(EMP(i,k)));
end;

ZZ = -L



