ARMA: Theory, Estimation, and Applications

• ARMA is a class of linear stochastic models that combine autoregressive and moving average components to model stationary processes.
• The model is expressed as a rational transfer function, facilitating forecasting, spectral analysis, and likelihood-based inference.
• Generalizations extend ARMA to spatial data, graph filters, and non-Gaussian settings, broadening its applicability in diverse scientific domains.

Searching arXiv for recent and foundational ARMA papers relevant to a comprehensive synthesis. ARMA, or autoregressive moving average, denotes a class of linear stochastic models in which a process is represented through an autoregressive polynomial and a moving-average polynomial, equivalently through a rational transfer function. In the standard stationary univariate setting, an ARMA$(p,q)$ model is written as

$$X_t - \sum_{i=1}^{p} \phi_i X_{t-i} = \varepsilon_t + \sum_{j=1}^{q} \theta_j \varepsilon_{t-j},$$

or, with the backshift operator $B$, as $\Phi(B)X_t=\Theta(B)\varepsilon_t$. Across the literature, the same rational-filter viewpoint extends beyond scalar time series to two-dimensional random fields, graph filters, quaternion-valued signals, and non-Gaussian or stochastic relaxations used for missing data, heavy tails, and biomedical imaging (Ganesh et al., 2024).

1. Classical formulation and spectral viewpoint

In its classical form, ARMA combines an autoregressive part, governed by coefficients $\phi_1,\dots,\phi_p$, with a moving-average part, governed by $\theta_1,\dots,\theta_q$. Using

$$\Phi(B)=1-\phi_1B-\cdots-\phi_pB^p,\qquad \Theta(B)=1+\theta_1B+\cdots+\theta_qB^q,$$

the model is $\Phi(B)X_t=\Theta(B)\varepsilon_t$, where $\varepsilon_t$ is white noise with $\mathbb{E}[\varepsilon_t]=0$ and variance $$X_t - \sum_{i=1}^{p} \phi_i X_{t-i} = \varepsilon_t + \sum_{j=1}^{q} \theta_j \varepsilon_{t-j},$$0. Stationarity requires the zeros of $$X_t - \sum_{i=1}^{p} \phi_i X_{t-i} = \varepsilon_t + \sum_{j=1}^{q} \theta_j \varepsilon_{t-j},$$1 to lie outside the closed unit disc, and invertibility requires the zeros of $$X_t - \sum_{i=1}^{p} \phi_i X_{t-i} = \varepsilon_t + \sum_{j=1}^{q} \theta_j \varepsilon_{t-j},$$2 to lie outside the unit disc (Ganesh et al., 2024).

This representation is equivalent to viewing ARMA as a linear time-invariant filter with rational transfer function

$$X_t - \sum_{i=1}^{p} \phi_i X_{t-i} = \varepsilon_t + \sum_{j=1}^{q} \theta_j \varepsilon_{t-j},$$3

In the moving-average representation $$X_t - \sum_{i=1}^{p} \phi_i X_{t-i} = \varepsilon_t + \sum_{j=1}^{q} \theta_j \varepsilon_{t-j},$$4, the generating function $$X_t - \sum_{i=1}^{p} \phi_i X_{t-i} = \varepsilon_t + \sum_{j=1}^{q} \theta_j \varepsilon_{t-j},$$5 determines the second-order structure, and the spectral density is

$$X_t - \sum_{i=1}^{p} \phi_i X_{t-i} = \varepsilon_t + \sum_{j=1}^{q} \theta_j \varepsilon_{t-j},$$6

For ARMA processes, $$X_t - \sum_{i=1}^{p} \phi_i X_{t-i} = \varepsilon_t + \sum_{j=1}^{q} \theta_j \varepsilon_{t-j},$$7 is rational; this is the core link between ARMA modeling and rational approximation on the unit circle (Ganesh et al., 2024).

From the forecasting perspective, the one-step-ahead predictor in the classical Gaussian formulation is

$$X_t - \sum_{i=1}^{p} \phi_i X_{t-i} = \varepsilon_t + \sum_{j=1}^{q} \theta_j \varepsilon_{t-j},$$8

with one-step predictive variance $$X_t - \sum_{i=1}^{p} \phi_i X_{t-i} = \varepsilon_t + \sum_{j=1}^{q} \theta_j \varepsilon_{t-j},$$9 when $B$0. This standard formulation underlies exact likelihood methods, state-space implementations, and later graphical-model and generalized extensions (Thiesson et al., 2012).

2. Likelihood, estimation, and inference

Classical ARMA estimation is typically performed by maximum likelihood, nonlinear optimization, or innovations/Kalman-filter methods. In widely used software, the exact Gaussian likelihood is commonly evaluated through a state-space representation and Kalman filtering, often initialized by conditional sum-of-squares (CSS). This workflow is effective but not benign: likelihood surfaces can be strongly multimodal, and single-start optimization can converge to suboptimal local maxima (Wheeler et al., 2023).

A systematic study of these optimization failures shows that standard single-start likelihood maximization improved in at least $B$1 of simulated cases when replaced by a root-based random-initialization algorithm, and in $B$2 of cases for $B$3 with $B$4. The same study reports that $B$5 of AIC tables built from single-start fits exhibit nested-model inconsistencies, and that profile confidence intervals provide superior confidence intervals to those based on the Fisher information matrix (Wheeler et al., 2023). This suggests that practical ARMA inference depends as much on optimization strategy and interval construction as on the likelihood itself.

A separate computational line replaces repeated exact ARMA likelihood evaluations by a likelihood-based AR approximation. The proposed approximation requires $B$6 work only once to precompute a small Champernowne matrix, after which repeated likelihood evaluations cost $B$7 in the series length for fixed approximation order. In most cases, the resulting estimates are identical to or very close to the exact maximum likelihood estimate, making high-level implementations in environments such as Mathematica, Matlab, and R practical for long series (McLeod et al., 2016).

3. Stochastic, graphical, and non-Gaussian generalizations

Expressing ARMA as a directed graphical model exposes a structural obstacle: in the classical model, $B$8 is deterministic given its parents, and this degeneracy makes EM effectively unusable. A stochastic relaxation replaces the deterministic observation equation by

$B$9

with small user-fixed $\Phi(B)X_t=\Theta(B)\varepsilon_t$0, yielding the $\Phi(B)X_t=\Theta(B)\varepsilon_t$1-ARMA and $\Phi(B)X_t=\Theta(B)\varepsilon_t$2-ARMA* models. This modification restores EM, supports missing-data inference, and accommodates cross predictors for multiple time series and nontemporal covariates (Thiesson et al., 2012).

Within this framework, EM alternates between Gaussian inference over latent innovations and regression-style M-step updates for $\Phi(B)X_t=\Theta(B)\varepsilon_t$3 and $\Phi(B)X_t=\Theta(B)\varepsilon_t$4. The graphical-model formulation also gives closed-form one-step predictive distributions of the form

$\Phi(B)X_t=\Theta(B)\varepsilon_t$5

with $\Phi(B)X_t=\Theta(B)\varepsilon_t$6. On complete economic time series, $\Phi(B)X_t=\Theta(B)\varepsilon_t$7-ARMA and $\Phi(B)X_t=\Theta(B)\varepsilon_t$8-ARMA* outperform classical ARMA in sequential predictive score; on US-Econo, for example, the average sequential predictive scores reported are $\Phi(B)X_t=\Theta(B)\varepsilon_t$9 for ARMA, $\phi_1,\dots,\phi_p$0 for smoothed ARMA, $\phi_1,\dots,\phi_p$1 for $\phi_1,\dots,\phi_p$2-ARMA, and $\phi_1,\dots,\phi_p$3 for $\phi_1,\dots,\phi_p$4-ARMA* (Thiesson et al., 2012).

When Gaussian innovations are inadequate, ARMA has also been extended to symmetric $\phi_1,\dots,\phi_p$5-stable noise and to ARMA–GARCH with stable innovations. In that setting, autocovariances are undefined for $\phi_1,\dots,\phi_p$6, so estimation is built on normalized autocovariation,

$\phi_1,\dots,\phi_p$7

a modified Yule–Walker step, and LAD regression in a modified Hannan–Rissanen procedure. For stable GARCH components, the paper proposes a modified empirical characteristic function method. The reported simulations show that LAD and the modified Hannan–Rissanen method substantially outperform least squares under heavy tails, and the methodology is applied to IBM log returns via an MA(1)–S$\phi_1,\dots,\phi_p$8S–GARCH(1,1) specification (Sathe et al., 2019).

4. Multidimensional and spatial ARMA

For images and spatial lattices, ARMA generalizes to two-dimensional random fields. In the breast-imaging formulation, an image is modeled as a field $\phi_1,\dots,\phi_p$9 on $\theta_1,\dots,\theta_q$0 with a total order

$\theta_1,\dots,\theta_q$1

and a 2D ARMA$\theta_1,\dots,\theta_q$2 model

$\theta_1,\dots,\theta_q$3

The model is estimated by a two-stage Yule–Walker least-squares procedure: first, a finite AR approximation to the inverse filter is estimated by 2D Yule–Walker; second, the resulting innovation field is inserted into a least-squares system for joint AR and MA estimation (0906.3722).

In that application, local ARMA parameters are used as texture features on $\theta_1,\dots,\theta_q$4 windows of $\theta_1,\dots,\theta_q$5 ultrasound images, and k-means with $\theta_1,\dots,\theta_q$6 segments healthy tissue, benign tumor, and cancerous tumor. The reported overall detection/classification accuracy is $\theta_1,\dots,\theta_q$7, with the method motivated by the claim that breast images can be accurately modeled by two-dimensional ARMA random fields and that the estimated coefficients provide a compact parametric descriptor of local tissue texture (0906.3722).

A different spatial generalization is the 2-D Rayleigh ARMA model for SAR amplitudes, where the observation at $\theta_1,\dots,\theta_q$8 is conditionally Rayleigh with mean $\theta_1,\dots,\theta_q$9, and a link-transformed mean $$\Phi(B)=1-\phi_1B-\cdots-\phi_pB^p,\qquad \Theta(B)=1+\theta_1B+\cdots+\theta_qB^q,$$0 follows a unilateral 2D ARMA recursion. The model is estimated by conditional maximum likelihood, with score recursion, Fisher information, and Wald tests derived explicitly. On CARABAS II forest data, the reported anomaly detector based on 2-D RARMA(1,1) detects $$\Phi(B)=1-\phi_1B-\cdots-\phi_pB^p,\qquad \Theta(B)=1+\theta_1B+\cdots+\theta_qB^q,$$1 of $$\Phi(B)=1-\phi_1B-\cdots-\phi_pB^p,\qquad \Theta(B)=1+\theta_1B+\cdots+\theta_qB^q,$$2 military vehicles with $$\Phi(B)=1-\phi_1B-\cdots-\phi_pB^p,\qquad \Theta(B)=1+\theta_1B+\cdots+\theta_qB^q,$$3 false alarms, whereas a Gaussian 2-D ARMA(1,1) detects $$\Phi(B)=1-\phi_1B-\cdots-\phi_pB^p,\qquad \Theta(B)=1+\theta_1B+\cdots+\theta_qB^q,$$4 vehicles with $$\Phi(B)=1-\phi_1B-\cdots-\phi_pB^p,\qquad \Theta(B)=1+\theta_1B+\cdots+\theta_qB^q,$$5 false alarms; for image modeling on the same data, the reported MSE and MAPE are $$\Phi(B)=1-\phi_1B-\cdots-\phi_pB^p,\qquad \Theta(B)=1+\theta_1B+\cdots+\theta_qB^q,$$6 and $$\Phi(B)=1-\phi_1B-\cdots-\phi_pB^p,\qquad \Theta(B)=1+\theta_1B+\cdots+\theta_qB^q,$$7 for RARMA versus $$\Phi(B)=1-\phi_1B-\cdots-\phi_pB^p,\qquad \Theta(B)=1+\theta_1B+\cdots+\theta_qB^q,$$8 and $$\Phi(B)=1-\phi_1B-\cdots-\phi_pB^p,\qquad \Theta(B)=1+\theta_1B+\cdots+\theta_qB^q,$$9 for Gaussian ARMA (Palm et al., 2022).

5. Graph, spectral, and neural ARMA filters

On graphs, ARMA becomes a rational filter in the eigenvalues of a graph operator. One formulation defines a first-order graph ARMA recursion

$\Phi(B)X_t=\Theta(B)\varepsilon_t$0

where $\Phi(B)X_t=\Theta(B)\varepsilon_t$1 is a shifted Laplacian, producing the graph-frequency response

$\Phi(B)X_t=\Theta(B)\varepsilon_t$2

Parallel and periodic ARMA$\Phi(B)X_t=\Theta(B)\varepsilon_t$3 constructions then realize higher-order rational graph filters. A central feature of this line of work is graph-independent design: coefficients are chosen for a spectral interval rather than for a particular graph, and the resulting filters are robust to changes in the signal and/or graph (Loukas et al., 2015).

A centralized design perspective writes the graph filter as

$\Phi(B)X_t=\Theta(B)\varepsilon_t$4

with frequency response

$\Phi(B)X_t=\Theta(B)\varepsilon_t$5

Because the true least-squares design problem is nonlinear, the literature develops Prony-inspired modified-error methods and iterative true-error minimization. On both synthetic and real graph data, the reported results show that ARMA graph filters outperform FIR graph filters in approximation accuracy and are suitable for interpolation, compression, and prediction (Liu et al., 2017).

In graph neural networks, ARMA filters appear as recursive convolutional layers. A Graph Convolutional Skip layer updates

$\Phi(B)X_t=\Theta(B)\varepsilon_t$6

and several such stacks approximate an ARMA$\Phi(B)X_t=\Theta(B)\varepsilon_t$7 spectral response. The reported experiments show improvements over polynomial graph filters on semi-supervised node classification, graph signal classification, graph classification, and graph regression, while preserving locality, trainability, and transferability across graphs (Bianchi et al., 2019). The same rational-filter idea has also been transplanted to dense prediction in Euclidean CNNs: ARMA layers replace ordinary convolutions by combining a moving-average convolution with output-to-output autoregressive couplings, yielding an adjustable receptive field governed by learnable AR coefficients and a stable re-parameterization for training (Su et al., 2020). More recently, ARMA graph convolutions have been used as the encoder in a contrastive and modularity-regularized biomedical node-classification framework, ARMA-C3, for ADNI, NIFD, BreastMNIST, PneumoniaMNIST, and liver ultrasound data (Abburi et al., 25 May 2026).

6. Online, distributed, and other operational settings

ARMA structure also governs distributed detection in temporally correlated environments. In the distributed detection setup where the signal is the impulse response of one ARMA filter and the noise is the output of another ARMA filter driven by white Gaussian noise, the paper extends the running consensus detector by whitening each local stream through fixed-order ARMA filters. The asymptotic theory yields two regimes: either local error probabilities decay exponentially fast to zero, or they converge to a strictly positive error floor. Necessary and sufficient conditions are given in terms of the poles and zeros of the composite ARMA models, and the threshold level influences the asymptotics in the positive-floor regime (Domingos et al., 2023).

Another extension replaces real or complex coefficients by quaternion-valued coefficients. The quaternion ARMA$\Phi(B)X_t=\Theta(B)\varepsilon_t$8 model

$\Phi(B)X_t=\Theta(B)\varepsilon_t$9

is learned online by reducing the hidden-noise problem to a full-information quaternion AR$\varepsilon_t$0 approximation. The resulting qARMA-QOGD and qARMA-QONS algorithms use quaternion GHR calculus and achieve logarithmic regret bounds; the reported synthetic experiments show that both approaches converge toward the optimal MSE, with the Newton-style method typically converging faster (Pu et al., 2019).

A common ambiguity is purely terminological: a separate system named “Arma” is a Byzantine Fault Tolerant consensus system that separates dissemination and validation of client transactions from consensus and orders only metadata of batches of transactions. Despite the name, that usage is unrelated to autoregressive moving-average modeling (Manevich et al., 2024).