# Julia functions to check for first or second order stochastic dominance
# between discrete random variables with finite support
# We define here two functions: depo(pa,xa,pb,xb), and deso(pa,xa,pb,xb)
# Ricard Torres <rtorres@itam.mx>
"""
function depo(pa,xa,pb,xb)\n
Check for first order stochastic dominance between the\n
finite discrete distributions (pa,xa) and (pb,xb)\n
pi are the probabilities attached to the different values\n
xi are values (in strict increasing order)\n
Output is: 0 if there is mutual domination (distributions coincide)\n
  1 if (pa,xa) dominates (pb,xb)\n 
 -1 if (pa,xa) is dominated by (pb,xb)\n
  2 if neither distribution dominates the other one.
"""
function depo(pa,xa,pb,xb)
   la = length(xa)
   lb = length(xb)
   if la != length(pa) || lb != length(pb)
      error("Values and probabilities must be vectors of the same length")
   end
   tol = 1e-10
   if max( abs(sum(pa)-1), abs(sum(pb)-1) ) > tol
      error("Probability vectors must add up to 1")
   end
   if xa != sort(xa) || xb != sort(xb) || xa != unique(xa) || xb != unique(xb)
      error("Vectors of values must be in strict ascending order")
   end
   xd = sort(unique(vcat(xa,xb))) # common support of xa and xb
   ld = length(xd)
   fd = zeros(ld+1) # prepend a zero
   ia = 1
   ib = 1
   for i = 2 : ld + 1
       fd[i] = fd[i-1] 
       if ia <= la && xd[i-1] == xa[ia]
          fd[i] += pa[ia]
          ia += 1
       end
       if ib <= lb && xd[i-1] == xb[ib]
          fd[i] -= pb[ib]
          ib += 1
       end
   end
   adominates = true
   bdominates = true
   for i = 2 : ld
       if fd[i] > tol
          adominates = false
       elseif fd[i] < - tol
          bdominates = false
       end
   end
   if adominates && bdominates 
      return 0   # mutual domination -> equality
   elseif adominates
      return 1   # (xa,pa) dominates (xb,pb)
   elseif bdominates
      return -1  # (xb,pb) dominates (xa,pa)
   else
      return 2   # neither one dominates the other
   end
end

"""
function deso(pa,xa,pb,xb)\n
Check for second order stochastic dominance between the\n
finite discrete distributions (pa,xa) and (pb,xb)\n
pi are the probabilities attached to the different values\n
xi are values (in strict increasing order)\n
Output is: 0 if there is mutual domination (distributions coincide)\n
  1 if (pa,xa) dominates (pb,xb)\n 
 -1 if (pa,xa) is dominated by (pb,xb)\n
  2 if neither distribution dominates the other one.
"""
function deso(pa,xa,pb,xb)
   la = length(xa)
   lb = length(xb)
   if la != length(pa) || lb != length(pb)
      error("Values and probabilities must be vectors of the same length")
   end
   tol = 1e-10
   if max( abs(sum(pa)-1), abs(sum(pb)-1) ) > tol
      error("Probability vectors must add up to 1")
   end
   if xa != sort(xa) || xb != sort(xb) || xa != unique(xa) || xb != unique(xb)
      error("Vectors of values must be in strict ascending order")
   end
   xd = sort(unique(vcat(xa,xb))) # common support of xa and xb
   ld = length(xd)
   fd = zeros(ld+1) # prepend a zero
   ia = 1
   ib = 1
   for i = 2 : ld + 1
       fd[i] = fd[i-1] 
       if ia <= la && xd[i-1] == xa[ia]
          fd[i] += pa[ia]
          ia += 1
       end
       if ib <= lb && xd[i-1] == xb[ib]
          fd[i] -= pb[ib]
          ib += 1
       end
   end
   adominates = true
   bdominates = true
   integd = zeros(ld)
   for i = 2 : ld
       integd[i] = integd[i-1] + fd[i] * ( xd[i] - xd[i-1] )
       if integd[i] > tol
          adominates = false
       elseif integd[i] < - tol
          bdominates = false
       end
   end
   if adominates && bdominates 
      return 0   # mutual domination -> equality
   elseif adominates
      return 1   # (xa,pa) dominates (xb,pb)
   elseif bdominates
      return -1  # (xb,pb) dominates (xa,pa)
   else
      return 2   # neither one dominates the other
   end
end