This paper introduces a rigorous class of two-dimensional Markov-switching autoregressive moving-average (2D MS-ARMA) models for spatial lattice data exhibiting regime-dependent dynamics. The switching mechanism is governed by a latent causal Markov random field that drives spatial transitions between regime-specific autoregressive and moving-average structures. We provide sufficient conditions for the existence of a strictly stationary solution through the top Lyapunov exponent associated with a sequence of random matrices obtained from a state-space representation constructed along the lexicographic order. For the first-order bidirectional specification, we derive explicit spectral conditions linking stationarity to the regime-dependent spectral radii. Sufficient conditions ensuring the existence of finite second-order moments are also provided. Parameter estimation is carried out using a variational expectation–maximization (VEM) algorithm based on a mean-field approximation of the posterior distribution of the hidden regimes. The E-step yields closed-form coordinate ascent updates, while the M-step relies on gradient-based numerical optimization with derivatives computed via recursive differentiation. Under increasing-domain asymptotics, we discuss the consistency and asymptotic behavior of the variational estimator. The proposed framework fills a methodological gap between classical one-dimensional Markov-switching ARMA models and spatial autoregressive structures by extending regime-switching theory to multi-indexed processes with rigorous probabilistic foundations. It provides a comprehensive basis for statistical inference, model diagnostics, and prediction in spatially heterogeneous environments.
Rashedi et al. (Wed,) studied this question.