High- and multi-dimensional array data are becoming increasingly available. They admit a natural representation as tensors and call for appropriate statistical tools. We propose a new linear autoregressive tensor process (ART) for tensor-valued data, that encompasses some well-known time series models as special cases. We study its properties and derive the associated impulse response function. We exploit the PARAFAC low-rank decomposition for providing a parsimonious parameterization and develop a Bayesian inference allowing for shrinking effects. We apply the ART model to time series of multilayer networks and study the propagation of shocks across nodes, layers and time.