Actually, co-skewness is represented by a rank 3 tensor, rather than a matrix.
I'm going to reproduce the formulation from Bhandari and Das, Options on portfolios with higher-order moments, but I'll add and omit some details.
The co-skewness tensor is
$$
S_{ijk} = E \left[ r_i \times r_j \times r_k \right] = \frac{1}{T} \sum_{t=1}^T r_i(t) \times r_j(t) \times r_k(t)
$$
where $r$ are asset returns over $T$ time periods.
Then, given portfolio weights $w$, mean asset returns $\mu$, covariance matrix $\Sigma$, and portfolio variance $\sigma_p^2(w) = w\prime \Sigma w$, we calculate moments:
$$\begin{eqnarray*}
m_1 & = & w\prime \mu \\
m_2 & = & \sigma_p^2 + m_1^2 \\
m_3 & = & \sum_{i=1}^N \sum_{j=1}^N \sum_{k=1}^N w_i w_j w_k S_{ijk}
\end{eqnarray*}$$
The portfolio skewness is then
$$
S_p = \frac{1}{\sigma_p^3} \left[ m_3 - 3m_2 m_1 + 2m_1^3 \right]
$$
In the case of a 6-asset portfolio, the co-skewness tensor will contain 216 components; however, due to symmetry, it only contains 56 independent components.
Therefore, it can be helpful to reformulate the portfolio skewness equation for computational efficiency. To do this, we can start with the definition of skewness for portfolio returns,
$$
S_p = \frac{1}{\sigma_p^3} E [ \left( \sum_{i=1}^N w_i r_i \right)^3 ] \quad ,
$$
and then apply the multinomial theorem to obtain the portfolio skewness in terms of only the independent components.
Update
- Especially for longer time series, the return moments should be
centered on the means, i.e., $r_i = R_i - \bar{R}_i$
- In the case of daily returns, $R_i(t) = \frac{P(t) -
P(t-1)}{P(t-1)}$, where $P(t)$ is the closing price at time $t$.
- Be sure the prices for returns are comparable from period to period.
For example, stock prices may need adjustments to account for
dividend payments. See this Q & A on return measurement for
more discussion.
note: I edited the equation for the co-skewness tensor above.