function es=es(x1,x2)

% This program carries out the Omnibus Test for Two Samples developed in Epps & Singleton 1986

% x1 and x2 are data for the two samples

x=[x1; x2];
T1=length(x1);
T2=length(x2);
T=length(x);

%These are the values for t1 and t2 recommended by the authors
t1=0.4;
t2=0.8;
t=[t1; t2];

%Here we choose the standardizing sigma for the t's according to p. 202

x_sorted=sort(x);
sigma_hat=0.5*((x_sorted(floor(T*0.25),:)+x_sorted(floor(T*0.25)+1,:))/2+(x_sorted(floor(T*0.75),:)+x_sorted(floor(T*0.75)-1,:))/2);
t_hat=t/sigma_hat;

%Creating the 2J x Ti matrices (J=2 and ti is the # obs for sample i), one for each of the two treatments, equation (3) in paper

J=length(t);
for i=1:T1
    for j=1:J
        g_x1m(1+(j-1)*J,i)=cos(t_hat(j)*x1(i));
        g_x1m(2+(j-1)*J,i)=sin(t_hat(j)*x1(i));
    end
end

for i=1:T2
    for j=1:J
        g_x2m(1+(j-1)*J,i)=cos(t_hat(j)*x2(i));
        g_x2m(2+(j-1)*J,i)=sin(t_hat(j)*x2(i));
    end
end

%Covariances S1_hat and S2_hat (the biased version is used: S1_hat_2 and S2_hat_2)

S1_hat=cov(g_x1m');
S2_hat=cov(g_x2m');

S1_hat_2=S1_hat*((T1-1)/T1);
S2_hat_2=S2_hat*((T2-1)/T2);

%Creating the gi's vectors

g1=(1/T1)*sum(g_x1m,2);
g2=(1/T2)*sum(g_x2m,2);

%Computing Omega hat

v1=T1/T;
v2=T2/T;

omega_hat=(1/v1)*S1_hat_2+(1/v2)*S2_hat_2;
omega_hat_inv=inv(omega_hat);

%Computing the Statistic
G2=(g1-g2);
if T1<25 & T2<25
    c_T1_T2=1/(1+(T1+T2)^(-0.45)+10.1*(T1^(-1.7)+T2^(-1.7)));
    W2=c_T1_T2*T*G2'*omega_hat_inv*G2;
else 
    W2=T*G2'*omega_hat_inv*G2;
end

%Computing p-value and critical values

r=rank(omega_hat);
p_value=1-chi2cdf(W2,r);
critical_value_10=chi2inv(0.90,r);
critical_value_5=chi2inv(0.95,r);
critical_value_1=chi2inv(0.99,r);

disp(['The value of the Epps-Singleton Statistic is: ' num2str(W2)])
disp(['The p-value is: ' num2str(p_value)])
disp(['The 10%, 5% and 1% critical values are: ' num2str(critical_value_10) ' ' num2str(critical_value_5) ' ' num2str(critical_value_1)])