Abstract: In several disciplines, it is common to find time series measured at irregular observational times. In particular, in astronomy there are a large number of surveys that gather information over irregular time gaps and in more than one passband. Some examples are Pan-STARRS, ZTF, and also the LSST. However, current commonly used time series models that estimate the time dependence in astronomical light curves consider the information of each band separately (e.g, CIAR, IAR, and CARMA models) disregarding the dependence that might exist between different passbands. In this paper, we propose a novel bivariate model for irregularly sampled time series, called the Bivariate Irregular Autoregressive (BIAR) model. The BIAR model assumes an autoregressive structure on each time series; it is stationary, and it allows to estimate the autocorrelation, the cross-correlation and the contemporary correlation between two unequally spaced time series. We implemented the BIAR model on light curves, in the g and r bands, obtained from the ZTF alerts processed by the ALeRCE broker. We show that if the light curves of the two bands are highly correlated, the model has more accurate forecast and prediction using the bivariate model than a similar method that uses only univariate information. Further, the estimated parameters of the BIAR are useful to characterize long-period variable stars and to distinguish between classes of stochastic objects, providing promising features that can be used for classification purposes.