%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%         Testing Stochastic Inequality given Independent Samples       %%
%%                             Matlab Code                               %%
%%                         disclaimer as usual                           %%
%%                          by Karl H. Schlag                            %%
%%                  August 4, 2008 - Copyrights reserved                 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%----------------------------------------------------------------------------------
%Reference: K. H. Schlag (2008, A New Method for Constructing Exact
%Tests without Making any Assumptions, Department of Economics and Business 
%Working Paper 1109, Universitat Pompeu Fabra
%----------------------------------------------------------------------------------

%Program for testing H0: P(Y2>Y1)-P(Y2<Y1)<= d vs H1: P(Y2>Y1)-P(Y2<Y1)> d
%given independent samples of Y1 and Y2
%need to have files test1 and test2 in workspace, 
%where test1 contains n1 observations of Y1 and test2 contains n2 observations of Y2

%to test H0: P(Y2>Y1)-P(Y2<Y1)>= d simply interchange Y1 and Y2 
%and test <=-d

close; clear; 

%enter parameters
T = 50000; %number of trials
d = 0; %cutoff for measure of stochastic inequality
theta = 0.2;
alpha = 0.05;
pseudoalpha = alpha*theta;

%read data
load('c:\karl\statistics\stochineq\datatest1');
load('c:\karl\statistics\stochineq\datatest2');

%calculating the sample mean (can be added to printout if desired)
n1 = length(test1);
mean1 = 0;
    for k=1:1:n1;
        mean1 = mean1 + test1(k,1)/n1;
    end;

n2 = length(test2);
mean2 = 0;
    for k=1:1:n2;
        mean2 = mean2 + test2(k,1)/n2;
    end;
    
%calculating unbiased estimate of P(Y2>Y1)-P(Y2<Y1)
sg = 0;
    for k=1:1:n1;
        d1=test2<test1(k);        
        d2=test2>test1(k);
        sg=sg+sum(d2)-sum(d1);
    end;
sg = sg /(n1*n2);
    
disp([' ']);    
disp(['estimate[P(Y_2>Y_1)-P(Y_2<Y_1)] = ', num2str(sg),'  alpha = ', num2str(alpha),'  theta = ', num2str(theta)]);

p=(1+d)/2; %threshold for binomial test
n=min(length(test1),length(test2));

%start loop
rj = 0;
for j=1:T;

%randomly pair data
    c1=test1(randperm(length(test1))'); %randomly permutate Y1
    c2=test2(randperm(length(test2))'); %randomly permutate Y2

%cut data so that we have an equal amount of data for Y1 and Y2
    c1=c1(1:n);
    c2=c2(1:n);

%adjust for d ne 0 by changing (0,0) and (1,1) sometimes 
    s1 = 0;
    s2 = 0;
    if d>0;
        for k=1:1:n;
            if c1(k) > c2(k);
                s1 = s1 + 1;
            elseif c1(k) < c2(k);
                s2 = s2 + 1;
            else
                q=rand;
                if q<d/(1+d);
                    s1 = s1 + 1;
                end;
            end;
        end;
    elseif d<0;
        for k=1:1:n;
            if c1(k) > c2(k);
                s1 = s1 + 1;
            elseif c1(k) < c2(k);
                s2 = s2 + 1;
            else
                q=rand;
                if q<-d/(1-d);
                    s2 = s2 + 1;
                end;
            end;
        end;
    else % d=0
        for k=1:1:n;
            if c1(k) > c2(k);
                s1 = s1 + 1;
            elseif c1(k) < c2(k);
                s2 = s2 + 1;
            end;
        end;
    end;
    
%Evaluation of randomized A test to transformed data (rejects if percentage
%of s2 among s1+s2 is significantly above p)
    h1 = 0;
    for k=s2:1:s1+s2;
        h1 = h1 + (p^k)*((1-p)^(s1+s2-k))*nchoosek(s1+s2,k);
    end;

    if h1 <= pseudoalpha; %(reject with probability 1)
        rj = rj + 1/T;
    else
        h = (p^s2)*((1-p)^s1)*nchoosek(s1+s2,s2);
        if h1 <= pseudoalpha + h; %(reject with positive probability)
            rj = rj + ((pseudoalpha - h1 + h) / h) /T ;
        end;
    end;

end; %of loop

if rj >= theta
    disp(['P(Y2>Y1)-P(Y2<Y1) <= ', num2str(d),' rejected  (rj=  ', num2str(rj),')']);    
else
    disp(['P(Y2>Y1)-P(Y2<Y1) <= ', num2str(d),' NOT rej.  (rj=  ', num2str(rj),')']);    
end;