obtaining 4th moment tensor under change of coordinates Suppose I have a random real-valued vector $x=x_1,\ldots,x_d$ and $M_{ijkl}=E[x_i x_j x_k x_l]$ and apply a change of coordinates $y=Ax$ where $A$ is an orthogonal matrix. How do I obtain $N_{ijkl}=E[y_i y_j y_k y_l]$?
 A: Using the Einstein summation convention, in which repeated indices are summed from $1$ through $d,$ and writing $A = (a_{ii^\prime}),$ we have
$$y_i = a_{ii^\prime}x_i^\prime.$$
Linearity of expectation yields
$$\begin{aligned}
N_{ijkl} &= E\left[y_iy_jy_ky_l\right] \\
&= E\left[a_{ii^\prime}x_i^\prime\, a_{jj^\prime}x_j^\prime\, a_{kk^\prime}x_k^\prime\, a_{ll^\prime}x_l^\prime\right]\\
&= a_{ii^\prime} a_{jj^\prime} a_{kk^\prime} a_{ll^\prime}E\left[x_i^\prime x_j^\prime x_k^\prime x_l^\prime\right] \\
&= a_{ii^\prime} a_{jj^\prime} a_{kk^\prime} a_{ll^\prime}M_{i^\prime j^\prime k^\prime l^\prime}.
\end{aligned}$$
This establishes $M$ as a (rank $4$) tensor.  Because there is no combination of the $a$'s (no two have common subscripts), no simplification is possible from the orthogonality assumption.
For more about algebraic manipulation of tensors, see e.g. https://en.wikipedia.org/wiki/Tensor#As_multidimensional_arrays or any advanced physics textbook.

Because this notation is not common among statisticians, and because some readers are uncomfortable with matrix formulas, here is a familiar example to relate these summations to a more familiar expression.  Consider the variance-covariance matrix, which (when all its components are well-defined and finite) is given by the formula
$$\Sigma_{ij}(\mathbf{x}) = E[x_ix_j] - E[x_i]E[x_j].$$
We seek the variance $\Sigma(\mathbf y)$ for the transformed random variable $\mathbf y=A\mathbf x$ for any $d\times d$ matrix $A = (a_{ij}).$
Again invoking linearity of expectation (which allows us to pull constant multiplicative factors like $a_{ii^\prime}$ out of the expectations), we may apply this series of simplifications:
$$\begin{aligned}
\Sigma(\mathbf{y})_{ij} &= E[y_iy_j] - E[y_i]E[y_j]\\
&= E[a_{ii^\prime}x_{i^\prime}\,a_{jj^\prime}x_{j^\prime}] - E[a_{ii^\prime}x_{i^\prime}]E[a_{jj^\prime}x_{j^\prime}]\\
&= a_{ii^\prime}a_{jj^\prime}E[x_{i^\prime}x_{j^\prime}] - a_{ii^\prime}E[x_{i^\prime}]\,a_{jj^\prime}E[x_{j^\prime}]\\
&= a_{ii^\prime}a_{jj^\prime}\left(E[x_{i^\prime}x_{j^\prime}] - E[x_{i^\prime}]E[x_{j^\prime}]\right)\\
&= a_{ii^\prime}\,a_{jj^\prime}\,\Sigma(\mathbf x)_{i^\prime j^\prime}\\
&= a_{ii^\prime}\,\Sigma(\mathbf x)_{i^\prime j^\prime}\,a^\prime_{j^\prime j}\\
&= (A\,\Sigma(\mathbf x)\,A^\prime)_{ij}.
\end{aligned}$$
The algebraic steps are simple, applying only the distributive and commutative laws of multiplication and then the definition of the transpose, $a^\prime_{j^\prime j} = a_{jj^\prime}.$  The last step applies the definition of matrix multiplication (namely, $(BC)_{ij}=b_{ik}c_{kj}$ for any compatible matrices $B=(b_{ik})$ and $C=(c_{kj})$) to express the $ij$ component of $\Sigma(\mathbf y)$ in terms of the matrices $A$ and $\Sigma(\mathbf x).$
Because two matrices are equal if and only if all their corresponding components are equal, we conclude

$$\Sigma(A\mathbf x) = A\,\Sigma(\mathbf x)\,A^\prime.$$

