Modeling stochastic variability in multiband time-series data

In preparation for the era of time-domain astronomy with upcoming large-scale surveys, we propose a state-space representation of a multivariate damped random walk process as a tool to analyze irregularly-spaced multifilter light curves with heteroscedastic measurement errors. We adopt a computationally efficient and scalable Kalman filtering approach to evaluate the likelihood function, leading to maximum O(k3n) complexity, where k is the number of available bands and n is the number of unique observation times across the k bands. This is a significant computational advantage over a commonly used univariate Gaussian process that can stack up all multiband light curves in one vector with maximum O(k3n3) complexity. Using such efficient likelihood computation, we provide both maximum likelihood estimates and Bayesian posterior samples of the model parameters. Three numerical illustrations are presented: (i) analyzing simulated five-band light curves for a comparison with independent singleband fits; (ii) analyzing five-band light curves of a quasar obtained from the Sloan Digital Sky Survey Stripe-82 to estimate short-term variability and timescale; (iii) analyzing gravitationally lensed g- and r-band light curves of Q0957+561 to infer the time delay. Two R packages, Rdrw and timedelay, are publicly available to fit the proposed models.