Computing the Co-Skewness matrix of portfolio

Dear all, I want to compute the skewness of portfolio. By reading loads of literature, I found that the formula of computing the skewness of portfolio is as following:

skewness_p=ω' M3 (ω⊗ω)

where M3=E[(r-μ) (r-μ)'⊗(r-μ)']={aijk}

the ⊗ denote the kronecker product, M3 is the co-skewness matrix, and r and u are the return matrix and average return matrix, respectively. However, I don't know how to decompose the expectation of the co-skewness matrix M3, thus I found a alternative way to calculate the M3, that's: The co-skewness matrix of M3 of dimensions (N,N2) can be represented by N Aijk matrixes (N,N) such that: M3=[ A1jk A1jk ... ANjk ] where j,k=1,…,N as well as for an index i. The individual elements of the matrix can be obtained as:

aijk=(1/T−1)∑i,j,k∑t [(Kt,i−μi)(Kt,j−μj)(Kt,k−μk)] .

the capital K denote the element in the TxN return matrix, T is the number of observations and N is the number of assets. According to the formula, I write down the following code:

M3=[];
for i=1:n
    S=[];
    for j=1:n
        for k=1:n
            u=0;
            for t=1:T 
                u=u+((P(t,i)-mean(P(:,i)))*(P(t,j)-mean(P(:,j))) ...
                    *(P(t,k)-mean(P(:,k))));
            end
            S(j,k)=u/(T-1); 
        end
    end
    M3=[M3 S];
end
sK=w*M3*kron(w',w')                % portfolio coskewness

P is my portfolio matrix.My problem comes out...There are 4 loops in this piece of code, and the speed is extremely low! I made a few simulation to approximate the time elapse. If I have 20 assets(n=20) in my portfolio, it will cost me 8 seconds, if I have 100 assets in my portfolio, it will cost me around 4 minutes. If I have 200 asset in my portfolio, it costs me 21 minutes! For the number of assets equal to 300, I run approximately 60 minutes and it still hasn't finished, so I have to abort it...The problem becomes unmanageable, I have 185,000 portfolios need to be evaluated, and around 17000 of them have more than 200 assets(n>=200). Within these 17000 portfolio, there are even more than 1200 portfolios with more than 1000 assets! This problem really makes me desperate! Does anyone know how can I improve this code? I reaaaaally reaaaally appreciate!!!

[Prach Arporniema]

I face the same problem with you. I can't solve it and I'm out of time now. For others who would process the higher moment optimization, I think matrix operations should minimize the calculation time. The loop makes it very slow. For myself, I only have to reduce the number of assets to fewer than 125.

Answer 1

[Peter Urbani]

For a Matlab implementation due to Pawel Lachowicz see Here

For an Excel and VBA implementation with open source code see Here

[Michele Gambera]

Lachowicz implementation is also not vectorized---same code as above. It even uses the same variable names.