function r = deso(pa, xa, pb, xb)
%% function r = deso(pa, xa, pb, xb)
%% Comprobar si hay dominancia estocastica de segundo orden
%% entre las distribuciones discretas (pa, xa) y (pb, xb)
%% pi son las probabilidades de cada valor
%% xi son los valores respectivos (en estricto orden ascendiente)
%% Resultado:
%% r =  0 si hay dominacion mutua (las distribuciones coinciden)
%% r =  1 si (pa, xa) domina a (pb, xb)
%% r = -1 si (pb, xb) domina a (pa, xa)
%% r =  2 si ninguna de las distribuciones domina a la otra
%%  Ricard Torres <rtorres@itam.mx>

   la = length(xa);
   lb = length(xb);
   if la ~= length(pa) || lb ~= length(pb)
      error('Valores y probabilidades deben tener la misma longitud');
   end
   tol = 1e-10; % tolerancia para comparar magnitudes
   if max( abs(sum(pa)-1), abs(sum(pb)-1) ) > tol
      error('Vectores de probabilidad deben sumar 1');
   end
   if any(xa ~= sort(unique(xa))) || any(xb ~= sort(unique(xb)))
      error('Vectores de valores deben estar en orden ascendiente estricto');
   end
   % soporte comun para xa y xb
   xd = sort(unique(horzcat(reshape(xa,1,la),reshape(xb,1,lb)))); 
   ld = length(xd);
   % Diferencia entre ambas distribuciones en el soporte comun
   fd = zeros(1,ld+1); % a~nadimos un cero al inicio
   ia = 1;
   ib = 1;
   for i = 2 : ld + 1
       fd(i) = fd(i-1); 
       if ia <= la && xd(i-1) == xa(ia)
          fd(i) = fd(i) + pa(ia);
          ia = ia + 1;
       end
       if ib <= lb && xd(i-1) == xb(ib)
          fd(i) = fd(i) - pb(ib);
          ib = ib + 1;
       end
   end
   adomina = true;
   bdomina = true;
   integd = zeros(1,ld); % integral de la diferencia
   for i = 2 : ld
       integd(i) = integd(i-1) + fd(i) * ( xd(i) - xd(i-1) );
       if integd(i) > tol
          adomina = false;
       elseif integd(i) < - tol
          bdomina = false;
       end
   end
   if adomina && bdomina 
      r = 0;   % dominacion mutua: igualdad
   elseif adomina
      r = 1;   % (pa,xa) domina a (pb,xb)
   elseif bdomina
      r = -1;  % (pb,xb) domina a (pa,xa)
   else
      r = 2;   % ninguna domina a la otra
   end
end