% This program replicates numerical example in Epps & Singleton 1986

%These are the entries for the two treatments
x1=[-0.15; 8.60; 5; 3.71; 4.29; 7.74; 2.48; 3.25; -1.15; 8.38];
x2=[2.55; 12.07; 0.46; 0.35; 2.69; -0.94; 1.73; 0.73; -0.35; -0.37];

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 (there is a mistake in the original paper, sigma_hat below
% is 2.05 and it is 1.95 in the paper

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);
sigma_hat=1.95;
t_hat=t/sigma_hat;

%Creating the 2J x 10 matrices (J=2 in this case), one for each of the two treatments, equation (3) in paper
J=2;

for i=1:length(x1)
    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:length(x2)
    for j=1:2
        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 (NOTE: THE S1 AND S2 REPORTED IN Epps and Singleton ARE NOT CORRECT...). However, these
%covariances S1_hat and S2_hat are not the ones used by Epps and Singleton. Matlab calculates the unbiased version of the variance
%and covariance by dividing through T-1 instead of T.

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

%Epps and Singleton RESULTS CAN BE REPRODUCED EXACTLY IF THE BIASED ESTIMATE OF VARIANCE AND COVARIANCE 
% (NOT NORMALIZED BY T-1) IS USED. To do this, we only rescale the covariances computed by Matlab accordingly:

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);

%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);

%G2
G2=(g1-g2);
W2=T*G2'*omega_hat_inv*G2;

%Small sample correction
c_T1_T2=1/(1+(T1+T2)^(-0.45)+10.1*(T1^(-1.7)+T2^(-1.7)));
W2_corrected=W2*c_T1_T2;

%Computing p-value and critical values
r=rank(omega_hat);
p_value=1-chi2cdf(W2_corrected,r);
critical_value_10=chi2inv(0.90,r);
critical_value_5=chi2inv(0.95,r);
critical_value_1=chi2inv(0.99,r);