The key here is in finding that $-$ for our application in finance $-$, the Kronecker product notation is 1) a way to shorten notation and 2) a function that is well represented in mathematical toolboxes such as
Assume there are $N=3$ random returns ($x_1,x_2,x_3$) and some weight vector $w$ of dimension $N\times 1$. The vector of expected returns with typical entry $E(x_i)$ has the dimension $\mu$ is $N\times1$, and $\mu^T\times w$ is a scalar ($1\times1$). The covariance matrix $\Sigma$ with typical entry $E((x_i-\mu_i)(x_j-\mu_j))$ has the dimension $N\times N$, and $w^T\times \Sigma \times w$ is again a scalar ($1\times1$). The co-skewness tensor $M_3$ with typical entry $E((x_i-\mu_i)(x_j-\mu_j)(x_k-\mu_k))$ has the dimension $N\times N\times N$, and again $w^TM_3\left(w\otimes w\right)$ is $1\times 1$. The same then holds for the kurtosis with $M_4$ of size $N\times N\times N \times N$.
Effectively, each statistics is simply an $K$-dimensional object, with $K=1,2,3,4$...
Now to the question: "How to compute the portfolio skewness?" (or kurtosis)
Assume your co-skewness matrix $M_3$ is stored as a three-dimensional array. You now right multiply the 3d-array ($N \times N \times N$) with a $N\times 1$ vector of portfolio weights $w$ $-$: the result is an $N\times N$ matrix! The next step is then the typical left/right multiplication with $w^T$ and $w$ et voilà: You have a scalar.
The algorithm (for the skewness) is:
- Build your $3$-dimensional skewness tensor $M_3$ with typical element as above
- For each entry $i$ in the last dimension, calculate $q_i=w^TM_3(.,.,i)w$. This yields a $N\times 1$ vector (by construction). Finally, $w^Tq$ gives you the portfolio skewness.
The same approach works for the portfolio kurtosis, i.e. reduce a 4d-array to 3d to 2d to 1d to scalar.