Is the third moment of an AR(1) dependent on $t$? Given an AR(1) process:
$$
X_t = \phi X_{t-1}+ \epsilon_t, \quad \epsilon\sim WN(0, \sigma^2)
$$
I know that if $|\phi|<1$, then the process is stationary (weakly). Thus, the first and second moment of $X_t$, $E(X_t)$ and $E(X_t^2)$, are constant. More specifically, they don't depend on $t$.
But, how about its third moment, $E(X_t^3)$, does it depend on $t$?

It may or may not be:
If $\epsilon_t$ is independent WN, the $MA(\infty)$ representation $X_t=\sum_{j=0}^\infty\phi^j\epsilon_{t-j}$ gives, for $|\phi|<1$,
$$
E(X_t^3)=\sum_{j=0}^\infty\phi^{3j}E(\epsilon_{t-j}^3),
$$
as pairs $\epsilon_i,\epsilon_j,\epsilon_k$ for which we do not have $i=j=k$ will yield terms of the form, e.g., $E(\epsilon_j^2)E(\epsilon_k)=0$.
If $E(\epsilon_{t-j}^3)$ is constant over time, and if we denote that quantity by $\gamma$, we obtain
$$
E(X_t^3)=\frac{\gamma}{1-\phi^3}
$$
A little illustration:
n <- 21000
k <- 2
epsilon <- rchisq(n, k)-k # a skewed mean zero distribution

phi <- 0.9
X <- arima.sim(model = list(ar=phi), n = n-1000, innov = epsilon[-(1:1000)], n.start = 1000, start.innov=epsilon[1:1000])

gamma <- k*(k+2)*(k+4) - 3*k^2*(k+2) + 3*k^3 - k^3 # 3rd moment of epsilon, see https://en.wikipedia.org/wiki/Chi-squared_distribution

gamma
mean(epsilon^3)

gamma/(1-phi^3)
mean(X^3)

