function [DH_stats DH_pvalue] = DoornikHansen(y)
% Conducts the Doornik-Hansen (2008) test
% INPUTS: PITs
% OUTPUTS: DH_stats - test statistics
%          DH_pvalue - the p-value associated with the test statistic
n = length(y);
for ll = 1:n
    y(ll)=min(y(ll),.999); % get rid of corner values
    y(ll)=max(y(ll),.001);
    z(ll)=norminv(y(ll));
end
[p n] = size(z);

% define moments
m2 = mean((z-mean(z)).^2);
m3 = mean((z-mean(z)).^3);
m4 = mean((z-mean(z)).^4);


sqrtb1 = m3.*(m2.^(-1.5));
b1 = sqrtb1.^2;
b2 = m4.*(m2.^(-2));

%Transformation of the skewness sqrtb1 by the D'Agostino procedure   
n2 = n.^2;
beta = 3*(n2 + 27*n - 70)*(n + 1)*(n + 3)/((n - 2)*(n + 5)*(n + 7)*(n + 9));
w2 = -1 + sqrt(2*(beta - 1));
d = 1/sqrt(log(sqrt(w2)));
y = sqrtb1*sqrt((w2 - 1)*(n + 1)*(n + 3)/(12*(n - 2)));
z1 = d*log(y + sqrt(y.^2 + 1)); %transformed skewness

%Transformation of kurtosiss from a gamma distribution to a chi-square
d = (n - 3)*(n + 1)*(n2 + 15*n - 4);
a = (n - 2)*(n + 5)*(n + 7)*(n2 + 27*n - 70)/(6*d);
c = (n - 7)*(n + 5)*(n + 7)*(n2 + 2*n - 5)/(6*d);
k = (n + 5)*(n + 7)*(n*n2 + 37*n2 + 11*n - 313) / (12*d);  
al = a + b1*c;
chi = (b2 - 1 - b1)*k*2;
z2 = (((chi./(2*al)).^(1/3)) - 1 + 1./(9*al)).*sqrt(9*al); %transformed kurtosis


DH_stats = z1*z1' + z2*z2'; %omnibus normality statistic
v = 2*p; %degrees of freedom
DH_pvalue = 1 - chi2cdf(DH_stats,v);