In this paper, we propose a vector autoregressive (VAR) model of order one for multivariate binary time series. Multivariate binary time series data are used in many fields such as biology and environmental sciences. However, modeling the dynamics in multiple binary time series is not an easy task. Most existing methods model the joint transition probabilities from marginals pairwisely for which the resulting cross-dependency may not be flexible enough. Our proposed model, Bernoulli VAR (BerVAR) model, is constructed using latent multivariate Bernoulli random vectors. The BerVAR model represents the instantaneous dependency between components via latent processes, and the autoregressive structure represents a switch between the hidden vectors depending on the past. We derive the mean and matrix-valued autocovariance functions for the BerVAR model analytically and propose a quasi-likelihood approach to estimate the model parameters. We prove that our estimator is consistent under mild conditions. We perform a simulation study to show the finite sample properties of the proposed estimators and to compare the prediction power with existing methods for binary time series. Finally, we fit our model to time series of drought events from different regions in Mexico to study the temporal dependence, in a given region and across different regions. By using the BerVAR model, we found that the cross-region dependence of drought events is stronger if a rain event preceded it.