function LM = LMAAstat(data)
% LM test based on Arellano-Valle and Azzalini(2008)
%   Detailed explanation goes here
N = size(data,1);
mu = mean(data);
Sigmainv = (cov(data,1))^(-1);
score = ((data - mu)*Sigmainv).^3;
I = eye(size(data,2));

% calculate the variance
for i = 1:size(data,2)
    for j = i:size(data,2)
        wij = I(i,:)* Sigmainv * I(:,j);
        wii = I(i,:)* Sigmainv * I(:,i);
        wjj = I(j,:)* Sigmainv * I(:,j);
        v_theta(i,j) = wij*(6*wij^2 + 9*wii*wjj);
        v_theta(j,i) = wij*(6*wij^2 + 9*wii*wjj);
    end
    v_thetamu(i,i) = 3*I(i,:)*Sigmainv*I(:,i);
end

v = v_theta - v_thetamu*Sigmainv*v_thetamu; % adjusted variance 

LM = N*mean(score)*v^-1*mean(score)';

end

