function p_value = hlntest(e1, e2, h)

if nargin < 2
    error('dmtest:TooFewInputs','At least two arguments are required');
end
if nargin < 3
   h = 1; 
end
if size(e1,1) ~= size(e2,1) || size(e1,2) ~= size(e2,2)
    error('dmtest:InvalidInput','Vectors should be of equal length');
end
if size(e1,2) > 1 || size(e2,2) > 1
    error('dmtest:InvalidInput','Input should have T rows and 1 column');
end
% Initialization
T = size(e1,1);
% Define the loss function from Harvey, Leybourne & Newbold (1998)
d = (e1-e2).*e1;
% Recalculate the variance of the loss differential, taking into account
% autocorrelation.
dMean = mean(d);
gamma0 = var(d);
if h > 1
    gamma = zeros(h-1,1);
    for i = 1:h-1
        sampleCov = cov(d(1+i:T),d(1:T-i),1);
        gamma(i) = sampleCov(2);
    end
    varD = gamma0 + 2*sum(gamma);
else
    varD = gamma0;
end
%k is calculated to adjust the statistic as per Harvey, Leybourne, and Newbold (1997) 
k = ((T+1-2*h+(((h)*(h-1))/T))/T)^(1/2);
% Retrieve the statistic HLN ~N(0,1)
    HLN = (dMean / sqrt ( (varD/T) ))*k;    
%P_VALUE is calculated    
    p_value = 2*tcdf(-abs(HLN),T-1);
end
    