My question is related to Moments of an AR(1) Process but with some subtle differences: I want to find the moments of more general ARMA(p,q) processes.
Context
I am working with a simple recurrence relation of the form:
\begin{equation} Y_i = aY_{i-1}+(1-a)X_i \quad i=1,2,\ldots n\tag{1} \end{equation}
where $Y_0=0$ and assuming $X_i$ are independent and identically distributed, but without assuming a distribution of the $X_i$. Given these conditions, I can then say:
\begin{equation} Y_n = (1-a)\sum_{i=1}^n a^{i-1}X_i \tag{2} \end{equation}
from which I want to identify the relationships between moments of $Y_n$ as a function of $X_n$which can be done by taking expectations of (2), to give:
\begin{align} \mathbf{E}\left[Y_n\right] &= \mathbf{E}\left[X_n\right]\\ \mathbf{Var}\left[Y_n\right] &= \frac{(1-a)^2}{1-a^2}\mathbf{Var}\left[X_n\right] = \frac{1-a}{1+a}\mathbf{Var}\left[X_n\right]\\ \mathbf{Skew}\left[Y_n\right] &= \frac{(1-a)^3}{1-a^3}\sqrt{\left(\frac{1-a^2}{(1-a)^2}\right)^3}\mathbf{Skew}\left[X_n\right] \end{align}
and so on, for $\mathbf{Kurt}\left[Y_n\right]$, $\mathbf{hyperskew}\left[Y_n\right]$ and $\mathbf{hyperkurt}\left[Y_n\right]$.
For a single recurrence relation such as (1) the algebra is tedious but do-able with care, but increases in complexity for the higher moments.
I am aware that the recurrence relation (1) is a form of AR model (an AR(1,0) perhaps, because the noise term is weighted and therefore there is a $\theta_0$ term). In ARMA-type formulation this is:
\begin{equation} (1-aB)Y_i = (1-a)X_i \tag{3} \end{equation}
This is a discretely-coincident form of the solution to the first-order differential equation:
\begin{equation} (1+aD)y(t) = x(t) \quad \mbox{for} \quad t>0 \end{equation}
where the inputs $x(t)$ are pulsed. Box and Jenkins describe this link between continuous and discrete-time models in their 1970 book and give a couple of examples in the Appendix to Chapter 4.
This allows us to recognise that $a=\mathrm{e}^{-1/k_a}$, which provides the appropriate averaging weights for the change from continuous to discrete time.
What I am looking for
If the output from (1) is then fed into a subsequent recurrence, and then the output from that into another, this can be written as an ARMA-type formulation. For example, if there are two recurrence relations (with parameters $a$ and $b$) then:
\begin{equation} \left(1-a B\right)\left(1-b B\right)y_m = \left(\frac{k_a}{k_a-k_b}\left(1-a\right)\left(1-bB\right)+\frac{k_b}{k_b-k_a}\left(1-b\right)\left(1-aB\right)\right)x_m \tag{4} \end{equation}
Is there a way to use this discretely-coincident form to find the moments? Or is there a way to express the moments of an ARMA(p,q) process such that I can use this discretely-coincident form to determine them?
The reason I am asking the question is because the moments for (1) are easily found from the discretely coincident form (3). However, this is more complex when the system is a cascade where the parameters ($a$, $b$) change between cascades, because of the weighting of terms that appears as shown in (4). It would be particularly useful to use just the discretely coincident form, because this also defines the covariance structure in the output.
This type of cascade from one linear relation to another is known in water resources as a cascade of linear reservoirs and is a core component of models that predict water movement (for resource planning, flooding, etc). I have found references that help me to analyse simple inputs-output relationships for a single linear system (such as that described in (1)), but not a cascade, and particularly not a cascade where the parameter $a$ changes along the cascade. I have made significant progress to derive the discretely-coincident form for such a cascade, but have not thus far been able to see how to identify the moment relationships as I have described. It looks to me as though this should be possible and I am wondering if there are some established methods for doing this given an ARMA formulation such as (4). I have searched for this through academic journals, this forum and elsewhere but have not found anything save the previous post referenced above.
Updated comment
I had thought that, given the discretely-coincident form in (4), this may be rearranged to:
\begin{equation} y_m = \left(\frac{k_a}{k_a-k_b}\frac{\left(1-a\right)}{\left(1-aB\right)}+\frac{k_b}{k_b-k_a}\frac{\left(1-b\right)}{\left(1-bB\right)}\right)x_m \tag{4} \end{equation}
and then this would allow the same logic as above to derive the moments, as the addition of those from AR(1) processes, but using the new weights. However, I have checked this numerically and it does not work.
However, I did compute some algebra using the recursion approach that assumed the $Y_i$ terms were then fed through a second recursion with parameter $b$ to produce a series $Z_i$. The variance from this can be shown to be:
\begin{equation} \mathbf{Var}\left[Z_n\right] = \mathbf{Var}\left[X_n\right] \left(\frac{1-a}{1+a}\right)\left(\frac{1-b}{1+b}\right)\left(\frac{1+ab}{1-ab}\right) \tag{5} \end{equation}
and I have determined this is correct through numerical simulations of a few different random variables fed through the process. But how would I get from (4) to (5) directly?
