  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.
