Mardia Skewness and Kurtosis for Multivariate Normal

Question

I have written a Matlab function to calculate Mardia's Skewness and Kurtosis. (S is the Covariance Matrix and X_bar is the mean vector.)

$$b_{1,p}=\frac{1}{n^{2}}\sum_{i,j=1}^{n}[(X_{i}-\bar{X})'S^{-1}(X_{j}-\bar{X})]^3$$

$$b_{2,p}=\frac{1}{n}\sum_{i=1}^{n}[(X_{i}-\bar{X})'S^{-1}(X_{i}-\bar{X})]^2$$

function [ b ] = Mardia(type, X )
%MARDIA Mardia's multivariate sample skewness and Kurtosis
%   Definition as presented in Schwager (2006)
%   LaTeX Def Skewness:
%       b_{1,p}=\frac{1}{n^{2}}\sum_{i,j=1}^{n}[(X_{i}-\bar{X})'S^{-1}(X_{i}-\bar{X})]^3
%   LaTeX Def Kurtosis:
%       b_{1,p}=\frac{1}{n^{2}}\sum_{i,j=1}^{n}[(X_{i}-\bar{X})'S^{-1}(X_{i}-\bar{X})]^4
X_bar = mean(X);
S = cov(X);
n = size(X,1);

% needs slow for loop due to memory restrictions
temp = 0;
switch type
    case 'Skewness'
        for i=1:n
            for j=1:n
                temp = temp + ( (X(i,:)-X_bar) * S^(-1) * (X(j,:)-X_bar)' )^3;
            end
            disp(strcat(num2str(i), '/', num2str(n)));
        end
        b = ( 1/(n^2)) * temp;
    case 'Kurtosis'
        for i=1:n
            temp = temp + ( (X(i,:)-X_bar) * S^(-1) * (X(i,:)-X_bar)' )^2;
        end
        b = ( 1/(n)) * temp;
end        
end

To check my calculation I have simulated a Multivariate Normal and calculated Mardia's Skewness (0.0077) and - Kurtosis (15.1989) (3 dimensional with 10000 draws).

What are Mardia's Multivariate Skewness and Kurtosis of a MVN?

I would also greatly appreciate any suggestions to speed up my code.