In this paper, we perform multivariate time series analysis from a Bayesian machine learning perspective through the proposed multivariate Bayesian time series (MBTS) model. The multivariate structure and the Bayesian framework allow the model to take advantage of the association structure among target series, select important features, and train the datadriven model at the same time. Extensive analyses on both simulated data and empirical data indicate that the MBTS model is able to, cover the true values of regression coefficients in 90% credible intervals, select the most important predictors, and boost the prediction accuracy with higher correlation in absolute value of the target series, and consistently yield superior performance over the univariate Bayesian structural time series (BSTS) model, the autoregressive integrated moving average with regression (ARIMAX) model, and the multivariate ARIMAX (MARIMAX) model, in one-step-ahead forecast and ten-steps-ahead forecast.