Multivariate Continuous-Time Autoregressive Moving Average Processes for Astronomical Multiband Time Series

Large-scale astronomical surveys provide unprecedented volumes of multivariate time-series observations obtained through multiple optical filters. We develop a structured multivariate continuous-time autoregressive moving average (MCARMA) framework for multi-band time series with irregular sampling, heteroscedastic measurement errors, and partially observed bands. The framework allows band-specific stochastic dynamics while modeling cross-band dependence through correlated Brownian driving processes, with state-space and spectral representations enabling likelihood-based inference and interpretation of the fitted stochastic dynamics. We develop a two-stage estimation procedure in which a numerically stabilized preliminary fit initializes subsequent maximum likelihood estimation. Simulations show that higher-order stochastic structure can be recovered when its characteristic features are adequately resolved, but can become weakly identifiable because of limited temporal resolution or near pole–zero cancellation. Joint multivariate estimation improves parameter recovery in 23 of 27 settings and spectral recovery in 26 of 27 settings relative to separate single-band fits. Three Sloan Digital Sky Survey Stripe 82 quasars, respectively favoring MCARMA(1,0), MCARMA(2,0), and MCARMA(2,1), illustrate how joint multiband modeling uses cross-band dependence to inform marginal dynamics and can yield different model-order and spectral inference. The methodology is implemented in the Python package mcarma.