function stat = Marginstat(data)

N = size(data,1);
mu = mean(data);
sigma = cov(data,1);

score = ((data-mu)./sqrt(diag(sigma))').^3;
omega = diag(diag(sigma))^(-1/2)*sigma*diag(diag(sigma))^(-1/2);

for i = 1:size(data,2)
    for j = i+1:size(data,2)
        v_theta(i,j) = 6*omega(i,j)^3 + 9*omega(i,j);
        v_theta(j,i) = v_theta(i,j);
    end
end
v_theta = v_theta + 15*eye(size(data,2));
v_adj = v_theta - 9*omega;
stat = N*mean(score)*v_adj^(-1)*mean(score)';

end

