Periodically Stationary GM Increments

• Periodically stationary generalized multiple increments is a framework that unifies cyclostationarity, multi-seasonality, fractional integration, and long-memory, enabling both integer and non-integer differencing over various lags.
• The methodology leverages spectral representations and Wiener–Kolmogorov filtering to transform complex multi-seasonal signals into a vector-stationary context for optimal estimation and interpolation under spectral uncertainty.
• It extends classical ARIMA-type models into robust minimax filtering applications, making it valuable for forecasting in environmental, economic, and communications data with overlapping or nested periodicities.

Periodically stationary generalized multiple increments (PS-GM-increments) form a unifying stochastic framework that merges cyclostationarity, multi-seasonality, (fractional) integration, and long-memory, applicable to both discrete-time and continuous-time settings. These structures generalize classical ARIMA/ARFIMA, SARFIMA, PARMA, and cyclostationary models by incorporating multiple seasonalities and allowing both integer and non-integer (fractional) differencing at potentially several lags, along with periodicities in mean and second moments. This formalism underpins minimax-robust estimation and interpolation of observed processes with multiple seasonalities and regimes in the presence of spectral uncertainty.

1. Formal Definition and Structure

Consider a discrete-time scalar sequence $\{\zeta(k)\}_{k\in\mathbb Z}$. For strictly positive integers $r$, seasonalities $\mathbf s=(s_1,\ldots,s_r)\in\mathbb N^r$, and orders $\mathbf d=(d_1,\ldots,d_r)\in\mathbb R^r$, define the generalized multiple (GM) difference operator: $$\Delta^{\mathbf d}(B)\equiv\prod_{j=1}^r(1-B^{s_j})^{d_j}$$ where $B$ is the backward-shift, $B\,\zeta(k)=\zeta(k-1)$. Each term $(1-B^{s_j})^{d_j}$ is expanded by the fractional binomial theorem, converging for $|B|=1$.

If $\zeta$ is $T$-periodic in second moments (i.e., cyclostationary), construct the $T$-dimensional vector $X(m)=(\zeta(mT+1),\ldots,\zeta(mT+T))^\top$. The increments $Y(m)=\Delta^{\mathbf{d}}(B)X(m)$ are PS-GM-increments if

$$\mathbb{E}[Y(m)]=\mathbb{E}[Y(m+1)], \quad \operatorname{Cov}(Y(m_1),Y(m_2))=\operatorname{Cov}(Y(m_1+1),Y(m_2+1))$$

for all $m_1$, $m_2$. The process is periodically stationary in its increments, and an equivalent vector-stationary process emerges in the “lifted” space.

In the continuous-time regime, let $X(t)$, $t\in\mathbb{R}$, and fix $T>0$; the $d$-th periodically correlated increment process is

$$\Delta X^{(d)}(t,T) = \sum_{\ell=0}^d (-1)^\ell \binom{d}{\ell} X(t-\ell T)$$

with the cyclostationarity property: $D^{(d)}(t+T,s+T) = D^{(d)}(t,s)$ for all $t$, $s$.

2. Spectral Representation and Covariance Structure

For discrete ($T$-vector) processes, the PS-GM-incremented sequence admits the following spectral representation (Luz et al., 2020, Luz et al., 10 Nov 2025, Luz et al., 2 Feb 2024): $$Y(m) = \int_{-\pi}^{\pi} e^{im\lambda} \Psi(e^{-i\lambda}) dZ(\lambda)$$ where $\Psi(e^{-i\lambda})$ encodes the transfer function arising from the multiplicative differencing,

$$\Psi(e^{-i\lambda})=\prod_{j=1}^r (1-e^{-is_j\lambda})^{-d_j}$$

and $dZ(\lambda)$ is a vector-valued orthogonal increment process.

The spectral density then factors as

$$f(\lambda) = \Psi(e^{-i\lambda}) \; f_0(\lambda) \; \Psi(e^{i\lambda})^*$$

where $f_0(\lambda)$ is a "short-memory" (innovation) spectral matrix, and each factor $(1-e^{-is_j\lambda})^{d_j}$ induces a fractional pole or zero at frequency $2\pi k/s_j$.

In scalar cases or under further factorization assumptions (Luz et al., 2023): $$f_\xi(\lambda) = |B^{d,P}(e^{-i\lambda})|^{-2}\, f(\lambda)$$ where $B^{d,P}$ is the composite multi-order, multi-seasonal difference operator.

For continuous-time periodically stationary increments, via expansion in an orthonormal basis of $L_2([0,T))$, define $$Y_{j,k}=\langle \Delta X^{(d)}(u+jT,T),e_k(u)\rangle$$ yielding an infinite-dimensional stationary sequence $Y_j = (Y_{j,1},Y_{j,2},\ldots)^\top$, and its spectral density matrix is

$$F(\lambda)_{k\ell} = (1-e^{-i\lambda T})^d \; f_{k\ell}(\lambda)\; \overline{(1-e^{-i\lambda T})^d}$$

where $f_{k\ell}(\lambda)$ is determined from the orthogonal increment decomposition (Luz et al., 2 Feb 2024).

3. Linear Forecasting, Filtering, and Interpolation

Given the observations $\zeta(k)+\eta(k)$ and specifying a linear functional $A(\zeta)=\sum_{k=0}^\infty a_k \zeta(k)$, the mean-square optimal linear estimator minimizing $\mathbb{E}|A(\zeta)-\hat{A}(\zeta)|^2$ is characterized in the frequency domain by the classical Wiener–Kolmogorov (Wiener–Hopf) formula. For a $T$-vector setting (Luz et al., 2020, Luz et al., 7 Nov 2025, Luz et al., 2 Feb 2024): $$h(\lambda) = \Psi(e^{-i\lambda}) \frac{\sum_k a_k e^{ik\lambda} f_0(\lambda)^* \Psi(e^{i\lambda})}{f(\lambda)}$$ subject to the constraint that $h$ is supported on the past spectral subspace $L_2((-\pi,0);\;f)$.

For noisy observations, the optimal filter's spectral characteristic becomes

$$H(e^{i\lambda}) = \frac{A(e^{-i\lambda}) f(\lambda)}{f(\lambda) + |B^{d,P}(e^{-i\lambda})|^2 g(\lambda)}$$

and the mean-square error is

$$\frac{1}{2\pi}\int_{-\pi}^{\pi} \frac{|A(e^{i\lambda})|^2 f(\lambda) |B^{d,P}(e^{-i\lambda})|^{-2} g(\lambda)}{f(\lambda) + |B^{d,P}(e^{-i\lambda})|^2 g(\lambda)} \, d\lambda$$

(Luz et al., 2023).

In the finite-interval ("interpolation") problem, block Toeplitz systems emerge, with explicit spectral formulae for mean-square error as matrix quadratic forms in the Fourier domain (Luz et al., 10 Nov 2025, Luz et al., 2021).

4. Minimax-Robust Filtering and Spectral Uncertainty

When the spectral densities $f(\lambda)$ (and/or $g(\lambda)$ for noise) are not known exactly, but belong to convex admissible sets (e.g., defined by energy, moment, trace, band-restriction, or $L_1$/trace neighborhoods), the estimation problem transitions to a minimax-robust setting (Luz et al., 2020, Luz et al., 7 Nov 2025, Luz et al., 2021, Luz et al., 2023, Luz et al., 2023): $$\min_h \max_{f\in \mathcal D} \mathbb{E}_{f} |A(\zeta) - \hat{A}_h(\zeta)|^2$$ or, for signal and noise jointly,

$$\min_H \max_{(f,g)\in \mathcal D_f\times\mathcal D_g}\mathbb{E}_{f,g}|A(\zeta) - H(\zeta+\eta)|^2$$

The least-favorable spectral densities $(f^*,g^*)$ that attain the maximal error for the optimal filter are characterized by systems of Lagrange-multiplier equations (KKT conditions) of the form (for almost every $\lambda$): $$|B(1-e^{-i\lambda T})^d|^2 [f^*(\lambda)+g^*(\lambda)]^{-2} = \Phi(\lambda)$$ with $\Phi(\lambda)$ depending on the constraint set and functional $A$.

Once $(f^*,g^*)$ are found, the minimax-robust filter substitutes these into the classical Wiener–Kolmogorov formula. Applications include explicit handling of spectral model misspecification, robust forecasting and interpolation for time series with uncertain periodic or seasonal structure, and design of communication filters tolerant to spectral uncertainty (Luz et al., 7 Nov 2025, Luz et al., 2020, Luz et al., 2021).

5. Continuous-Time and Higher-Order Generalization

The theory extends seamlessly to continuous time $X(t)$, with increments of the form $(1-B_T)^d X(t)$ and periodic structure in mean and covariance: $R_{\Delta X}(t+T, s+T) = R_{\Delta X}(t,s)$ This periodic correlation allows reduction to an infinite-dimensional vector-valued stationary sequence, enabling explicit spectral representations and optimal filter synthesis as in the discrete case (Luz et al., 2 Feb 2024, Luz et al., 2023).

Generalized multiple increments, allowing for $r$ lags and orders, result in increment sequences of the form

$$\chi_{\overline{\mu},\overline{s}}^{(d)}(\xi(m)) = \sum_{l_1=0}^{d_1}\cdots\sum_{l_r=0}^{d_r}(-1)^{l_1+\dots+l_r} \binom{d_1}{l_1} \cdots \binom{d_r}{l_r} \xi(m - \sum_{j=1}^r \mu_j s_j l_j)$$

with cyclostationary second-order and spectral structure, explicitly encoded in block matrix form.

This generalization is essential for analyzing signals or time series with overlapping or nested periodicities (e.g., multiple seasonal cycles in environmental, economic, or climatological data), and for handling higher-order differencing necessary for stationarization or spectral pole placement (Luz et al., 10 Nov 2025, Luz et al., 7 Nov 2025).

6. Implementation and Unified Modelling Context

Implementation of PS-GM-increment forecasting or filtering for a given dataset proceeds by:

1. Model identification: Estimating $r$, $(s_1,\dots,s_r)$ (seasonal periods), and orders $(d_1,\dots,d_r)$, possibly via frequency-domain or time-domain diagnostics.
2. Spectral estimation: Computing or prescribing $f(\lambda)$, $g(\lambda)$ via periodogram or model-based methods.
3. Increment expansion: Expressing increments in the required vectorized ($T$-block) form.
4. Solving block Toeplitz systems: Calculating filter coefficients via spectral integrals or solving matrix Wiener–Hopf equations, as dictated by boundary conditions (filtering, prediction, interpolation).
5. Minimax-robustification (if required): Formulating and solving the Lagrange/KKT system yielding least-favorable $(f^*,g^*)$ within prescribed admissible classes.
6. Filter synthesis: Constructing the filter in the frequency domain and, by inverse FFT or spectral integration, forming time-domain coefficients for application to observed increments (Luz et al., 7 Nov 2025, Luz et al., 2021).

This framework includes, as special cases, SARIMA, SARFIMA, PSARIMA, Gegenbauer/Cyclical ARFIMA, multivariate periodic autoregressions, and models for cointegrated sequences.

7. Connections, Limitations, and Applications

The PS-GM-increment methodology not only embodies classical cyclostationarity [Hurd & Miamee] and multi-seasonal long memory [Dudek et al.], but also establishes a general structure for minimax-robust signal processing and time series estimation under both certainty and uncertainty of spectral properties.

Limitations involve the requirement of square integrability, the necessity for invertible spectral densities in some derivations, and the challenge of high computational complexity for large $T$ or many seasonalities.

Applications of this theory include robust long-term forecasting of economic, environmental, and climatological time series exhibiting multi-periodic and long-range dependence, filtering of periodic communications signals, and sensitivity analysis to parametric spectral mis-specification. The approach is well-suited for scenarios where classical stationary or even cyclostationary models are insufficient to capture multiple interacting periodicities or long-memory/seasonal regime changes (Luz et al., 10 Nov 2025, Luz et al., 7 Nov 2025, Luz et al., 2023, Luz et al., 2 Feb 2024).