Consider a sequence of IID random variables $X_1,X_2,X_3,...$ from a common distribution with mean $\mu$, variance $\sigma^2$, skewness $\gamma$ and kurtosis $\kappa$ (all finite). The mean and variance of the sample mean are well-known:
$$\mathbb{E}(\bar{X}_n) = \mu \quad \quad \quad \quad \quad \mathbb{V}(\bar{X}_n) = \frac{\sigma^2}{n}.$$
What are the corresponding formulae for the skewness and kurtosis of the sample mean?
$$\begin{align} \boxed{ \quad \quad \ \ \ \mathbb{E}(\bar{X}_n) = \mu \quad \quad \quad \quad \quad \quad \quad \ \ \mathbb{V}(\bar{X}_n) = \frac{\sigma^2}{n}, \\[12pt] \quad \mathbb{Skew}(\bar{X}_n) = \frac{\gamma}{\sqrt{n}} \quad \quad \quad \quad \quad \mathbb{Kurt}(\bar{X}_n) = 3 + \frac{\kappa - 3}{n}. \quad \\} \end{align}$$
The mean, variance, skewness and kurtosis of the sample mean are shown in the box above. These formulae are valid for any case where the underlying values are IID with finite kurtosis. It is simple to confirm that $\mathbb{Skew}(\bar{X}_n) \rightarrow 0$ and $\mathbb{Kurt}(\bar{X}_n) \rightarrow 3$ as $n \rightarrow \infty$, which means that the sample mean is asymptotically unskewed and mesokurtic. This is also implied by the classical central limit theorem, which ensures that the standardised sample mean converges in distribution to the normal distribution.
Proof via cumulants: The simplest way to prove these results is via moment/cumulant generating functions. Suppose we let $m_X$ and $K_X$ denote the moment generating function and cumulant generating function for the underlying IID random variables in the sequence. It is simple to show that $m_{\bar{X}_n}(t) = m_X(t/n)^n$ so the cumulant generating function for the sample mean has the form:
$$K_{\bar{X}_n}(t) = n K_X(t/n).$$
Consequently, the cumulants of the sample mean are related to the cumulants of the underlying random variables by:
$$\begin{align} \bar{\kappa}_r \equiv \frac{d^r}{dt^r} K_{\bar{X}_n} (t) \Bigg|_{t=0} &= n \frac{d^r }{dt^r} K_X(t/n) \Bigg|_{t=0} \\[6pt] &= \frac{1}{n^{r-1}} K_X^{(r)}(t/n) \Bigg|_{t=0} \\[6pt] &= \frac{1}{n^{r-1}} K_X^{(r)}(0) \\[6pt] &= \frac{\kappa_r}{n^{r-1}}. \\[6pt] \end{align}$$
Using the relationship of the cumulants to the moments of interest, we then have:
$$\begin{align} \mathbb{V}(\bar{X}_n) &= \bar{\kappa}_2 \\[12pt] &= \frac{\kappa_2}{n} \\[6pt] &= \frac{\sigma^2}{n}, \\[6pt] \mathbb{Skew}(\bar{X}_n) &= \frac{\bar{\kappa}_3}{\bar{\kappa}_2^{3/2}} \\[6pt] &= \frac{\kappa_3 / n^2}{(\kappa_2/n)^{3/2}} \\[6pt] &= \frac{1}{\sqrt{n}} \cdot \frac{\kappa_3}{\kappa_2^{3/2}} \\[6pt] &= \frac{\gamma}{\sqrt{n}}, \\[6pt] \mathbb{Kurt}(\bar{X}_n) &= \frac{\bar{\kappa}_4 + 3 \bar{\kappa}_2^2}{\bar{\kappa}_2^2} \\[6pt] &= \frac{\kappa_4/n^3 + 3 (\kappa_2/n)^2}{(\kappa_2/n)^2} \\[6pt] &= \frac{\kappa_4/n + 3 \kappa_2^2}{\kappa_2^2} \\[6pt] &= \frac{(\kappa \sigma^4 - 3\sigma^4)/n + 3 \sigma^4}{\sigma^4} \\[6pt] &= 3 + \frac{\kappa - 3}{n}. \\[6pt] \end{align}$$
This method gives the results of interest, and it can also be generalise to give corresponding results for higher-order moments. As can be seen, the higher-order results are particularly simple, but the corresponding relationships for higher-order moments get messy once you get to high order.
Proof via expansion to raw moments: An alternative method of deriving these results is to expand the relevant central moments for the sample mean and simplify down to raw moments of the underlying random variables. Let $Y_i \equiv X_i - \mu$ and note that these random variables have mean zero, but have the same higher-order moments as $X_i$. The relevant higher-order central moments for the sample mean are:$^\dagger$
$$\begin{align} \mathbb{E}((\bar{X}_n - \mu)^3) &= \mathbb{E} \Bigg( \bigg( \frac{1}{n} \sum_{i=1}^n (X_i - \mu) \bigg)^3 \Bigg) \\[6pt] &= \mathbb{E} \Bigg( \bigg( \frac{1}{n} \sum_{i=1}^n Y_i \bigg)^3 \Bigg) \\[6pt] &= \frac{1}{n^3} \cdot \mathbb{E} \Bigg( \sum_{i=1}^n \sum_{j=1}^n \sum_{k=1}^n Y_i Y_j Y_k \Bigg) \\[8pt] &= \frac{1}{n^3} \cdot \mathbb{E} \Bigg( \sum_{i} Y_i^3 + \sum_{i \neq j} Y_i^2 Y_j + \sum_{i \neq j \neq k} Y_i Y_j Y_k \Bigg) \\[8pt] &= \frac{1}{n^3} \cdot \sum_{i} \mathbb{E}(Y_i^3) \\[12pt] &= \frac{1}{n^3} \cdot \sum_{i} \gamma \sigma^3 \\[12pt] &= \frac{1}{n^3} \cdot n \gamma \sigma^3 \\[12pt] &= \frac{\gamma}{n^2} \cdot \sigma^3, \\[12pt] \mathbb{E}((\bar{X}_n - \mu)^4) &= \mathbb{E} \Bigg( \bigg( \frac{1}{n} \sum_{i=1}^n (X_i - \mu) \bigg)^4 \Bigg) \\[6pt] &= \mathbb{E} \Bigg( \bigg( \frac{1}{n} \sum_{i=1}^n Y_i \bigg)^4 \Bigg) \\[6pt] &= \frac{1}{n^4} \cdot \mathbb{E} \Bigg( \sum_{i=1}^n \sum_{j=1}^n \sum_{k=1}^n \sum_{l=1}^n Y_i Y_j Y_k Y_l \Bigg) \\[6pt] &= \frac{1}{n^4} \cdot \mathbb{E} \Bigg( \sum_{i} Y_i^4 + \sum_{i \neq j} Y_i^3 Y_j + \sum_{i \neq j} Y_i^2 Y_j^2 \\[12pt] &\quad \quad \quad \quad \quad \quad + \sum_{i \neq j \neq k} Y_i^2 Y_j Y_k + \sum_{i \neq j \neq k \neq l} Y_i Y_j Y_k Y_l \Bigg) \\[6pt] &= \frac{1}{n^4} \cdot \Bigg[ \sum_{i} \mathbb{E}(Y_i^4) + \sum_{i \neq j} \mathbb{E}(Y_i^2) \mathbb{E}(Y_j^2) \Bigg] \\[6pt] &= \frac{1}{n^4} \cdot \Bigg[ \sum_{i} (\kappa \sigma^4) + \sum_{i \neq j} \sigma^4 \Bigg] \\[6pt] &= \frac{1}{n^4} \cdot \Bigg[ n \kappa \sigma^4 + 3n(n-1) \sigma^4 \Bigg] \\[6pt] &= \frac{(\kappa + 3(n-1)) \sigma^4}{n^3} \\[6pt] &= \frac{3n + (\kappa - 3)}{n^3} \cdot \sigma^4. \\[6pt] \end{align}$$
Consequently, the skewness and kurtosis of the sample mean are given respectively by:
$$\begin{align} \mathbb{Skew}(\bar{X}_n) &= \frac{\mathbb{E}((\bar{X}_n - \mu)^3)}{\mathbb{V}(\bar{X}_n)^{3/2}} \\[6pt] &= \frac{\gamma}{n^2} \cdot \sigma^3 \bigg/ \frac{\sigma^3}{n^{3/2}} \\[6pt] &= \frac{\gamma}{\sqrt{n}}, \\[12pt] \mathbb{Kurt}(\bar{X}_n) &= \frac{\mathbb{E}((\bar{X}_n - \mu)^4)}{\mathbb{V}(\bar{X}_n)^2} \\[6pt] &= \frac{3n + (\kappa - 3)}{n^3} \cdot \sigma^4 \bigg/\frac{\sigma^4}{n^2} \quad \quad \quad \quad \quad \quad \quad \\[6pt] &= \frac{3n + (\kappa - 3)}{n} \\[12pt] &= 3 + \frac{\kappa - 3}{n}. \\[12pt] \end{align}$$
$^\dagger$ We have used a slight abuse of notation in the range of the summations; when we refer to e.g., $i \neq j \neq k$ we use this as shorthand for the set of indices where all the indices are distinct. That is, we do not read these inequalities in their strict meaning, but read them as if the inequality were intended to be transitive.
The calculation amounts to removing the linear term in the cumulant generating function (log characteristic function) of the distribution and simply replacing its argument $t$ by the factor $t\sqrt{n}/\sigma$ needed to standardize the sum, afterwards multiplying everything by $n$ to account for the $n$ iid variables comprising the sum.
Let $\phi$ be the characteristic function of that common distribution with mean $\mu,$ standard deviation $\sigma,$ skewness $\gamma,$ and finite kurtosis $\kappa.$ Then because
$$\log \phi(t) = i \mu (t\sigma) - \frac{1}{2} (t\sigma)^2 - \frac{i\gamma}{6} (t\sigma)^3 + \frac{\kappa - 3}{24}(t\sigma)^4 + o(t^4),$$
the log characteristic function of the standardized sample mean $Z = \sqrt{n}(\bar X - \mu)/\sigma$ is
$$\begin{aligned} \log \phi_Z(t) &= n \left(\log \phi\left(\frac{t\sqrt{n}}{\sigma}\right) - i\mu\left(\frac{t\sqrt{n}}{\sigma}\right)\right)\\ &= -\frac{1}{2}\left(\frac{t}{\sqrt{n}}\right)^2 - n\frac{i\gamma}{6} \left(\frac{t}{\sqrt{n}}\right)^3 + n\frac{\kappa - 3}{24}\left(\frac{t}{\sqrt{n}}\right)^4 + o(t^4).\\ \end{aligned}$$
Comparing like powers of $t$ shows
$$\gamma_Z=\gamma/\sqrt{n},\ \kappa_Z - 3 = (\kappa-3)/n.$$
The first equation is standard--it's usually taken as the definition of skewness and kurtosis--while the second is a direct (and simple) consequence of the independence of the $X_i$ in the sample together with the rules for how $\phi$ transforms under recentering and rescaling. The rest of this post provides all details for any readers who might be unacquainted with this use of characteristic functions.
The characteristic function of a random variable $X$ is defined to be
$$\phi_X: \mathbb{R}\to \mathbb{C},\ \phi_X(t) = E\left[e^{itX}\right].$$
It always exists because $E\left[\big|e^{itX}\big|\right] \le E[1]=1$ demonstrates absolute convergence of the integrand $e^{itX}.$
Characteristic functions are useful for understanding (positive integral) moments $\mu^\prime(k) = E[X^k]$ because, when $E[X^{k}]$ exists and is finite, an application of Taylor's Theorem to the exponential shows
$$\begin{aligned} \phi_X(t) &= E\left[(1 + itX + (itX)^2/2! + \cdots + (itX)^k/k! + o(t^{k+1})\right]\\ &= 1 + \frac{i\mu_X(1)^\prime}{1}t^1 - \frac{\mu_X(2)^\prime}{2}t^2 - \frac{\mu_X(3)^\prime}{6}t^3 + \cdots + \frac{i^k \mu_X(k)^\prime}{k!}t^k + o(t^{k}). \end{aligned}$$
Thus, $\phi$ has a Maclaurin expansion $\phi_X(t) = a_0 + a_1/1\, t^1 + a_2/2!\, t^2 + \cdots + a_k/k!\, t^k$ (partial power series at $0$) and the moments of $X$ can be read directly from the coefficients $a_j:$ $\mu_X(j)^\prime = a_j(i)^j.$
A sample is defined to be a collection of $n$ independent random variables $X_i,$ $i=1,2,\ldots,n$ having a common distribution, whence all the $\phi_{X_i}$ are the same function $\phi.$ The sample mean is (also by definition)
$$\bar X = \left(X_1 + X_2 + \cdots + X_n\right)/n = X_1/n + X_2/n + \cdots + X_n/n.$$
Therefore, because the exponential of a sum of numbers is the product of their exponentials and expectations of products of independent random variables are the products of their expectations,
$$\phi_{\bar X}(t) = E\left[e^{it\bar X}\right] = E\left[e^{itX_1/n}\right]\, E\left[e^{itX_2/n}\right]\cdots E\left[e^{itX_1/n}\right] = \phi\left(\frac{t}{n}\right)^n.$$
The central moments of $X,$ where $\mu = \mu_X^\prime(1)$ exists and is finite, are defined as
$$\mu_X(k) = E\left[(X-\mu)^k\right] = \mu^\prime_{X-\mu}(k).$$
Consequently they may be found from the characteristic function of $X-\mu,$ which can be related to that of $X$ via
$$\phi_{X-\mu}(t) = E\left[e^{it(X-\mu)}\right] = E\left[e^{-it\mu}\,e^{itX}\right] = e^{-it\mu}\phi_X(t).$$
The relationship between the characteristic functions of $X$ and $X/\sigma,$ for any positive number $\sigma,$ is obtained directly from the definitions as
$$\phi_{X/\sigma}(t) = E\left[e^{itX/\sigma}\right] = \phi_X\left(\frac{t}{\sigma}\right).$$
The power relation between $\phi_{\bar X}$ and $\phi$ suggests working with the logarithms of these functions, because that would make this a direct proportion,
$$\psi_{\bar X}(t) = \log \phi_{\bar X}(t) = \log\left[\phi\left(\frac{t}{n}\right)^n\right] = n \log \phi\left(\frac{t}{n}\right) = n\psi(t).$$
Since, for $|s|\lt 1$ the Taylor series for $\log(1+s) = s -s^2/2 + s^3/3 - \cdots$ converges absolutely, we easily obtain the series
$$\psi_{X-\mu}(t) = \log\left(1 + \left[ - \frac{\mu_X(2)}{2}t^2 - \frac{\mu_X(3)}{6}t^3 + \frac{\mu_X(4)}{24}t^4 + o(t^4)\right]\right)$$
by setting $s = -\mu(2)(X)t^2 + \cdots + o(t^4)$ and computing
$$\begin{aligned}s^2/2&= \frac{1}{2}(-\mu_X(2)/2)^2 t^2 + o(t^4);\\ s^k/k! &= o(t^4)\end{aligned}$$
for all $k \gt 2.$
This gives
$$\psi_{X-\mu}(t) = 1 + s - s^2/2 + o(t^4) = - \frac{\mu_X(2)}{2}t^2 - i\frac{\mu_X(3)}{6}t^3 + \frac{\mu_X(4) - 3\mu_X(2)^2}{24}t^4 + o(t^4).$$
Upon scaling this by $1/\sigma = \mu_X(2)^{-1/2}$ it simplifies to
$$\psi_{(X-\mu)/\sigma}(t) = \log \phi_{(X-\mu)/\sigma}(t) =-\frac{1}{2}t^2 - \frac{i\gamma_X}{6}t^3 + \frac{\kappa_X - 3}{24}t^4 + o(t^4)$$
where $\gamma_X$ is, by definition, the skewness of $X$ and $\kappa_X$ is its kurtosis.
What are the first four central moments of the sample mean?
Computer algebra systems are particularly adept at this sort of algebra munching. Here I am using the mathStatica package for Mathematica where $s_r = \sum _{i=1}^n X^r$. The first say 7 central moments of the sample mean $\frac{s_1}{n}$ are:
where:
$\mu_r$ denotes the $r^\text{th}$ central moment of the population.
The calculation takes just .1 of a second for all of them to be calculated live.
