Let $\mathbf{Y}$ be a random vector. Are $k$th moments of $\mathbf{Y}$ considered? I am self-learning on linear model theory right now, and one thing I find surprising is that although $\mathbb{E}[\mathbf{Y}]$ is defined for a random vector $\mathbf{Y} = \begin{bmatrix}
y_1 \\
y_2 \\
\vdots \\
y_n\end{bmatrix}$, there is no mention of further moments besides the covariance matrix.
Google searching hasn't turned up much. Are $k$th (raw) moments of $\mathbf{Y}$ considered, or is there a different idea I don't know about?
I am learning from the text Plane Answers to Complex Questions (the TOC starts in p. 17 of the linked file). By "considered," what I mean is is there such a thing as $\mathbb{E}\left[\mathbf{Y}^k\right]$, and if so, how would such a concept be defined? The book I have only covers the first raw moment, and I find it a bit strange that there is no mention of how to define $\mathbb{E}\left[\mathbf{Y}^k\right]$ given my experience in univariate probability, nor do I have the expertise to define it.
Furthermore, if $\mathbb{E}\left[\mathbf{Y}^k\right]$ isn't defined, is there perhaps a related concept that I don't know about that is used instead?
 A: The proper analog of univariate moments in a multivariate setting is to view the exponent $\mathbf{k} = (k_1, k_2, \ldots, k_n)$ as a vector, too.  The exponential notation with vector bases and vector exponents is a shorthand for the product,
$$\mathbf{y}^\mathbf{k} = y_1^{k_1} y_2 ^{k_2} \cdots y_n^{k_n}.$$
For any such vector $\mathbf{k}$, the (raw) $\mathbf{k}^\text{th}$ moment of the random variable $\mathbf{Y}$ is defined to be
$$\mu_\mathbf{k} = \mathbb{E}\left(\mathbf{Y}^\mathbf{k}\right).$$
To motivate such a definition, consider a univariate moment of a linear function of $\mathbf{Y}$:
$$\mathbb{E}\left(\left(\lambda_1 Y_1 + \cdots + \lambda_n Y_n\right)^m\right) = \sum_\mathbf{k} \binom{m}{\mathbf{k}}\lambda^\mathbf{k} \mu_\mathbf{k}$$
where the sum occurs over all $\mathbf{k}$ whose components are whole non-negative numbers summing to $m$ and $\binom{m}{\mathbf{k}} = m!/(k_1!k_2!\cdots k_n!)$ are the multinomial coefficients.  The appearance of the multivariate moments on the right hand side shows why they are natural and important generalizations of the univariate moments.
These show up all the time.  For instance, the covariance between $Y_i$ and $Y_j$ is none other than
$$\text{Cov}(Y_i, Y_j)  = 
\mathbb{E}(Y_i Y_j)- \mathbb{E}(Y_i)\mathbb{E}(Y_j) = \mu_{\mathbf{k}_i + \mathbf{k}_j} - \mu_{\mathbf{k}_i}\mu_{\mathbf{k}_j}$$
where $\mathbf{k}_i$ and $\mathbf{k}_j$ are the indicator vectors with zeros in all but one place and a one in the indicated location.  (The same formula elegantly yields the variance of $Y_i$ when $i=j$.)
There are natural generalizations of all univariate moment concepts to the multivariate setting: a moment generating function, cumulants, a cumulant generating function, central moments, a characteristic function, and algebraic and analytical relationships among them all.
Reference
Alan Stuart and J. Keith Ord, Kendall's Advanced Theory of Statistics, Fifth Edition.  Oxford University Press, 1987: Volume I, Chapter 3, Moments and Cumulants.
A: In addition to @whuber's points 
1) I am not sure what linear model theory entails but remember that in linear models we are generally dealing with normal random variables which have 0 skew and 0 kurtosis.
2) More generally, the question is of the form "How precise is precise?". If I want to describe IID samples I could say I only want the mean. Alternatively I could say I want the mean and the errors in the means. An even more detailed alternative would be means, errors in the means and errors in the errors in the means. From this pattern you can see how the higher moments keep mounting. There is no real solution to this problem so people generally stop at level 2 (i.e., mean and variance). That is not to say the higher moments are useless. In fact, for problems involving fat-tailed distributions these problems become relevant
