  Multivariate spatio-temporal data arise more and more frequently in a wide
range of applications; however, there are relatively few general statistical
methods that can readily use that incorporate spatial, temporal and variable
dependencies simultaneously. In this paper, we propose a new approach to
represent non-parametrically the linear dependence structure of a multivariate
spatio-temporal process in terms of latent common factors. The matrix structure
of observations from the multivariate spatio-temporal process is well reserved
through the matrix factor model configuration. The spatial loading functions
are estimated non-parametrically by sieve approximation and the variable
loading matrix is estimated via an eigen-analysis of a symmetric non-negative
definite matrix. Though factor decomposition along the space mode is similar to
the low-rank approximation methods in spatial statistics, the fundamental
difference is that the low-dimensional structure is completely unknown in our
setting. Additionally, our method accommodates non-stationarity over space. The
estimated loading functions facilitate spatial prediction. For temporal
forecasting, we preserve the matrix structure of observations at each time
point by utilizing the matrix autoregressive model of order one MAR(1).
Asymptotic properties of the proposed methods are established. Performance of
the proposed method is investigated on both synthetic and real datasets
