The autoregressive analysis of stationary multivariate time series includes construction of a time domain model in the form of a multivariate stochastic difference equation and its transformation into a spectral matrix to study frequency domain properties. The time series is treated as a linear system with one output process and one or several inputs. This approach provides an explicit time domain description of the time series as a function of the past of output and input components. In the bivariate case, the frequency domain results include spectral densities, coherent spectrum, coherence function, and gain and phase factors. The coherence function shows the degree of linear interdependence, the coherent spectrum presents the contribution of the input component to the output, the gain factor consists of amplification coefficients, and the phase factor describes lags between the time series. All these quantities are frequency dependent. The Granger causality concept and feedbacks within the system are discussed as a part of bivariate time and frequency domain analysis. Some comments are given regarding the software for parametric analysis of multivariate time series. The use of cross-correlation coefficient and regression equation for describing dependence between time series or for reconstructing them is shown to be incorrect.