%% Funcion de Matlab/Octave para comprobar dominacion de Lorenz
%% entre distribuciones finitas de ingresos/riqueza
%% Ricard Torres <rtorres@itam.mx>
function c = lorenz_domin(na,ya,nb,yb)
% function c = lorenz_domin(na,ya,nb,yb)
% Comprobar si hay dominacion en el sentido de Lorenz entre 
% distribuciones finitas de riqueza (es suficiente comprobarlo 
% para los puntos que definen las curvas).
% Retorno: 
%  0: curvas coinciden 
%  1: (na,ya) domina a (nb,yb)
% -1: (nb,yb) domina a (na,ya) 
%  2: no hay dominacion
   la = length(na);
   lb = length(nb);
   if la ~= length(ya) || lb ~= length(yb)
      error('Vectores de poblacion e ingreso deben tener la misma longitud');
   end
   [Fa,Sa] = lorenz_disc(na,ya);
   [Fb,Sb] = lorenz_disc(nb,yb);
   F = sort(unique(horzcat(Fa(2:end-1),Fb(2:end-1))));
   lF = length(F);
   lora = interp1(Fa,Sa,F);
   lorb = interp1(Fb,Sb,F);
   tol = 1e-15; % Tolerancia para comparaciones
   a_encima = false;
   b_encima = false;
   for i = 1 : lF
       if ( lora(i) > lorb(i) && lora(i)-lorb(i) > tol )
          a_encima = true;
       elseif ( lorb(i) > lora(i) && lorb(i)-lora(i) > tol  )
          b_encima = true;
       end
   end
   if ( a_encima )
       if ( b_encima )
          c = 2;
          %disp('No hay dominacion de Lorenz')
       else
          c = 1;
          %disp('Primera serie domina a la segunda')
       end
   else
       if ( b_encima )
          c = -1;
          %disp('Segunda serie domina a la primera')
       else
          c = 0;
          %disp('Ambas curvas de Lorenz coinciden')
       end
   end
end

% Ejemplo:
% ya = [10 20 40 50 80];
% na = [40 40 80 20 20];
% yb = ya;
% nb = [30 40 90 30 10];
% c = lorenz_domin(na,ya,nb,yb)
% [Fa,Sa] = lorenz_disc(na,ya);
% [Fb,Sb] = lorenz_disc(nb,yb);
% plot(Fa,Sa,'-or',Fb,Sb,'-ob',[0,1],[0,1],'--g');
% xlabel('Frecuencias cumulativas poblacion'); 
% ylabel('Proporciones cumulativas ingreso');
% title('Dominacion de Lorenz');
% legend('Curva Lorenz Serie A','Curva Lorenz Series B','Diagonal','Location','NorthWest');