If an ARMA model formulation is written in infinite moving-average form: \begin{equation} X_t = C\left(B\right)\epsilon_t \quad \mbox{with} \quad C\left(B\right)=C_0+C_1B+C_2B^2 + \ldots \end{equation} where $\epsilon_t$ is zero-mean white noise with variance $\sigma_{\epsilon}^2$ but with non-zero skewness and moment coefficient of kurtosis $\neq3$ then this paper suggests that if we define $C_0=1$ then:
\begin{equation} \frac{\mathrm{E}\left[X_t^3\right]}{\left(\mathrm{E}\left[X_t^2\right]\right)^\frac{3}{2}}= \frac{\displaystyle\sum_{i=0}^{\infty}C_i^3}{\left(\displaystyle\sum_{i=0}^{\infty}C_i^2\right)^{\frac{3}{2}}}\mathbf{skew}\left[\epsilon\right] \quad \mbox{and} \quad \frac{\mathrm{E}\left[X_t^4\right]}{\left(\mathrm{E}\left[X_t^2\right]\right)^2}= \frac{\displaystyle\sum_{i=0}^{\infty}C_i^4}{\left(\displaystyle\sum_{i=0}^{\infty}C_i^2\right)^2}\mathbf{kurt}\left[\epsilon\right] + \frac{6\displaystyle\sum_{i=0}^{\infty} \sum_{j=i+1}^{\infty}C_i^2C_j^2}{\left(\displaystyle\sum_{i=0}^{\infty}C_i^2\right)^2} \end{equation}
It is not clear to me how these expressions were derived in two respects - first, the appearance of the additional term in the kurtosis formulation suggests some binomial expansion where single powers vanish, but I can't figure out what that expansion would be; and, (2) it is also not clear how the summation limits are determined in the kurtosis formulation. Any assistance/guidance gratefully received.
I assume that I need to expand a bracket that is equivalent to $\mathbf{E}\left[\left(X_t-\mu_{X_t}\right)^n\right]$ where $n=3$ for the skewness and $n=4$ for the kurtosis - and so on. I have tried doing this for $X_t = \mu + \left(1+\sum_{i=1}^{\infty}C_i \epsilon_{t-i}\right)$, but I am still unable to replicate the results above.
Update: I can generate very similar results to those given above if I start with an $\mathrm{AR}(1)$ process written in recursive form: \begin{equation} Y_i = aY_{i-1} + (1-a)X_i \end{equation} and successively substitute to gain an infinite moving average process. Expanding powers of this recursive expression allows me to obtain expressions for the higher moments, and taking expectations then allows me to find the expressions for the moments in terms of the parameter $a$.
The variance of the output $Y_n$ can be found from the second central moment of the above recurrence relation as: \begin{align}\nonumber \left(Y_i-\mu_{Y_i}\right)^2 &=\bigg(a\left(Y_{i-1}-\mu_{Y_{i-1}}\right)+\left(1-a\right)\left(X_i-\mu_{X_i}\right)\bigg)^2\\\label{eqn:lr_recurrence_variance_expansion} &=a^2\left(Y_{i-1}-\mu_{Y_{i-1}}\right)^2 + 2a\left(1-a\right)\left(Y_{i-1}-\mu_{Y_{i-1}}\right)\left(X_i-\mu_{X_i}\right) + \left(1-a\right)^2\left(X_i-\mu_{X_i}\right)^2 \end{align} Taking expectations gives: \begin{align}\nonumber \mathbf{E}\left[\left(Y_i-\mu_{Y_i}\right)^2\right] = a^2\mathbf{E}\left[\left(Y_{i-1}-\mu_{Y_{i-1}}\right)^2\right] &+ 2a\left(1-a\right)\mathbf{E}\left[\left(Y_{i-1}-\mu_{Y_{i-1}}\right)\right]\mathbf{E}\left[\left(X_i-\mu_{X_i}\right)\right]\\ &+ \left(1-a\right)^2\mathbf{E}\left[\left(X_i-\mu_{X_i}\right)^2\right] \end{align} and then note that: \begin{equation}\label{eqn:lr_moments_firstexpansionterms} \mathbf{E}\left[\left(Y_{i-1}-\mu_{Y_{i-1}}\right)\right]=\mathbf{E}\left[\left(X_i-\mu_{X_i}\right)\right]=0 \end{equation} such that: \begin{equation} \mathbf{E}\left[\left(Y_i-\mu_{Y_i}\right)^2\right] = a^2\mathbf{E}\left[\left(Y_{i-1}-\mu_{Y_{i-1}}\right)^2\right] + \left(1-a\right)^2\mathbf{E}\left[\left(X_i-\mu_{X_i}\right)^2\right] \end{equation} This is another recursion equation, for which we can determine the solution as: \begin{equation} \mathbf{E}\left[\left(Y_n-\mu_{Y_n}\right)^2\right] = a^{2n}\mathbf{E}\left[\left(Y_0-\mu_{Y_0}\right)^2\right] + \left(1-a\right)^2\left\{\sum_{i=0}^{n-1} a^{2i}\right\}\mathbf{E}\left[\left(X_n-\mu_{X_n}\right)^2\right] \end{equation} Given that $Y_0=0$ we can say that: \begin{align}\nonumber \mathbf{var}\left[Y_n\right] &= \left(1-a\right)^2\left\{\sum_{i=0}^{n-1} a^{2i}\right\}\mathbf{var}\left[X_n\right]\\\nonumber &= \left(\frac{\left(1-a\right)^2\left(1-a^{2n}\right)}{1-a^2}\right)\mathbf{var}\left[X_n\right]\\\label{eqn:lr_variance_output} &= \frac{\left(1-a\right)^2}{1-a^2}\mathbf{var}\left[X_n\right] \quad \mbox{or} \quad \left(\frac{1-a}{1+a}\right)\mathbf{var}\left[X_n\right] \quad \mbox{as} \quad n\rightarrow \infty \end{align}
In a similar way to the second central moment (variance) for the skewness we start from the third central moment \begin{align}\nonumber \left(Y_i-\mu_{Y_i}\right)^3 =a^3\left(Y_{i-1}-\mu_{Y_{i-1}}\right)^3 &+ 3a^2\left(1-a\right)\left(Y_{i-1}-\mu_{Y_{i-1}}\right)^2\left(X_i-\mu_{X_i}\right)\\\nonumber &+3a\left(1-a\right)^2\left(Y_{i-1}-\mu_{Y_{i-1}}\right)\left(X_i-\mu_{X_i}\right)^2\\\label{eqn:lr_recurrence_3rdcmoment_expansion} &+\left(1-a\right)^3\left(X_i-\mu_{X_i}\right)^3 \end{align} then take expectations and remove terms that equate to $0$ to give: \begin{equation} \mathbf{E}\left[\left(Y_i-\mu_{Y_i}\right)^3\right] = a^3\mathbf{E}\left[\left(Y_{i-1}-\mu_{Y_{i-1}}\right)^3\right] + \left(1-a\right)^3\mathbf{E}\left[\left(X_i-\mu_{X_i}\right)^3\right] \end{equation} This is, once again, a recursion equation for which the solution can be determined as: \begin{equation} \mathbf{E}\left[\left(Y_n-\mu_{Y_n}\right)^3\right] = a^{3n}\mathbf{E}\left[\left(Y_0-\mu_{Y_0}\right)^3\right] + \left(1-a\right)^3\left\{\sum_{i=0}^{n-1} a^{3i}\right\}\mathbf{E}\left[\left(X_n-\mu_{X_n}\right)^3\right] \end{equation} Once again, given that $Y_0=0$ we can say that: \begin{align}\nonumber \mathbf{E}\left[\left(Y_n-\mu_{Y_n}\right)^3\right] &= \left(1-a\right)^3\left\{\sum_{i=0}^{n-1} a^{3i}\right\}\mathbf{E}\left[\left(X_n-\mu_{X_n}\right)^3\right]\\ &=\left(\frac{\left(1-a\right)^3\left(1-a^{3n}\right)}{1-a^3}\right)\mathbf{E}\left[\left(X_n-\mu_{X_n}\right)^3\right]\\\label{eqn:lr_3rdcmoment} &= \frac{\left(1-a\right)^3}{1-a^3} \mathbf{E}\left[\left(X_n-\mu_{X_n}\right)^3\right] \quad \mbox{as} \quad n\rightarrow \infty \end{align} which gives the third central moment, from which we can determine the skewness also using the result for the variance to give: \begin{align}\nonumber \mathbf{skew}\left[Y_n\right] &= \frac{\mathbf{E}\left[\left(Y_n-\mu_{Y_n}\right)^3\right]}{\left(\mathbf{E}\left[\left(Y_n-\mu_{Y_n}\right)^2\right]\right)^\frac{3}{2}}\\ &=\left\{\frac{\left(1-a\right)^3}{1-a^{3}}\sqrt{\left(\frac{1+a}{1-a}\right)^3}\right\}\mathbf{skew}\left[X_n\right] \end{align}
The same process can be used to determine the kurtosis as:
\begin{align} \mathbf{kurt}\left[Y_n\right] &= \frac{\displaystyle 6a^2\left(1-a\right)^4\left\{\sum_{i=0}^{n-1} a^{4i}\right\}\left\{\sum_{i=0}^{n-1} a^{2i}\right\}\left(\mathbf{E}\left[\left(X_n-\mu_{X_n}\right)^2\right]\right)^2}{\displaystyle \left(1-a\right)^4\left\{\sum_{i=0}^{n-1} a^{2i}\right\}^2\left(\mathbf{E}\left[\left(X_n-\mu_{X_n}\right)^2\right]\right)^2}\\ &+\frac{\displaystyle \left(1-a\right)^4\left\{\sum_{i=0}^{n-1} a^{4i}\right\}\mathbf{E}\left[\left(X_n-\mu_{X_n}\right)^4\right]}{\displaystyle \left(1-a\right)^4\left\{\sum_{i=0}^{n-1} a^{2i}\right\}^2\left(\mathbf{E}\left[\left(X_n-\mu_{X_n}\right)^2\right]\right)^2} \end{align}
which becomes: \begin{align}\nonumber \mathbf{kurt}\left[Y_n\right] &= 6a^2 \left(\frac{1-a^{4n}}{1-a^4}\right)\left(\frac{1-a^2}{1-a^{2n}}\right) + \left(\frac{1-a^{4n}}{1-a^4}\right)\left(\frac{1-a^2}{1-a^{2n}}\right)^2\mathbf{kurt}\left[X_n\right]\\ &= \frac{6a^2}{1+a^2} + \left(\frac{1-a^2}{1+a^2}\right)\mathbf{kurt}\left[X_n\right] \quad \mbox{as} \quad n\rightarrow \infty \end{align}
It is clear that there is some correspondence between the two sets of results, but I can't figure out how the results from a specific infinite moving average formulation of an $\mathrm{AR}(1)$ process can be mapped to the more general $\mathrm{ARMA}$ solution proposed above.
Any thoughts or pointers welcome. N.B. I added the “cumulants” tag because I realised I might use them to solve the problem - but doing so gives different answers to these published solutions!
