I am working with complicated sets of time series, where an observed system output is the sum of a number of individual ARMA, error and noise series. I am therefore looking at ways to identify these individual series (or at least to speculate on the likely series from various characteristics of the aggregate autocovariance and spectral density functions).
As part of this, I need to use the autocovariance generating function as an easy way to convert infinite moving-average representations of ARMA models into their respective autocovariance and power spectral density functions.
However, I am struggling to understand the logic behind the autocovariance generating function. I have read and noted this post (How to retrieve information of autocovariance from autocovariance generating function?) but it doesn't explain why the derivation can be made. I have similarly read Hamilton's "Time Series" and followed his derivation of the ACVGF - but it is still not clear why.
The most complete derivation I could find went something like this:
For example, suppose $X_t$ is a linear process such that it can be written: \begin{equation} X_t = \sum_{i=0}^\infty \psi_iW_{t-i} = \psi\left(B\right)W_t \end{equation} and then the autocovariance function: \begin{equation} \gamma_h=\mathrm{Cov}\left(X_t,\,X_{t+h}\right)=\mathrm{E}\left[\sum_{i=0}^\infty \psi_iW_{t-i}\sum_{j=0}^\infty \psi_jW_{t+h-j}\right]=\sigma_w^2\sum_{i=0}^\infty \psi_i\psi_{i+h} \end{equation} Define the autocovariance generating function as: \begin{equation} g_{\gamma}\left(B\right)=\sum_{h=-\infty}^\infty \gamma_h B^h \end{equation} then: \begin{equation} g_{\gamma}\left(B\right)=\sigma_w^2\sum_{h=-\infty}^\infty\sum_{i=0}^\infty \psi_i\psi_{i+h} B^h=\sigma_w^2\sum_{i=0}^\infty\sum_{j=0}^\infty \psi_i\psi_{j}B^{j-i}=\sigma_w^2\sum_{i=0}^\infty \psi_iB^{-i}\sum_{j=0}^\infty \psi_j B^j=\sigma_w^2\psi\left(B^{-1}\right)\psi\left(B\right) \end{equation} Then if we note that: \begin{equation} g_{\gamma}\left(B\right)=\sum_{h=-\infty}^\infty \gamma_hB^h \end{equation} and: \begin{equation} f\left(\omega\right) = \sum_{h=-\infty}^\infty \gamma_h \mathrm{e}^{-2\pi i \omega h} = \gamma\left(\mathrm{e}^{-2\pi i \omega}\right)=\sigma_w^2\psi\left(\mathrm{e}^{-2\pi i \omega}\right)\psi\left(\mathrm{e}^{2\pi i \omega}\right)=\sigma_w^2\left|\psi\left(\mathrm{e}^{2\pi i \omega}\right)\right|^2 \end{equation} thus, if the form of the linear process can be expressed as an infinite moving average, then the spectrum can be automatically determined.
How can we just define the autocovariance generating function as that? What allows us to use the backshift operator (B) interchangeably with a complex valued scalar? How does this allow the spectrum to be directly determined by this substitution?
I assume that there are some implicit steps to the logic that are missing from this derivation - what are they and how can they be explained in straightforward terms?
