Modeling Stochastic Variability in Multi-Band Time Series Data

In preparation for the era of the 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 multi-filter 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(k^3n)$ 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 multi-band light curves in one vector with maximum $O(k^3n^3)$ 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 single-band fits; (ii) analyzing five-band light curves of a quasar obtained from the Sloan Digital Sky Survey (SDSS) Stripe~82 to estimate the 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.