Spatial and temporal correlations which affect the signal measured in functional MRI (fMRI) are usually not considered simultaneously (i.e., as non-independent random processes) in statistical methods dedicated to detecting cerebral activation. We propose a new method for modeling the covariance of a stationary spatio-temporal random process and apply this approach to fMRI data analysis. For doing so, we introduce a multivariate regression model which takes simultaneously the spatial and temporal correlations into account. We show that an experimental variogram of the regression error process can be fitted to a valid nonseparable spatio-temporal covariance model. This yields a more robust estimation of the intrinsic spatio-temporal covariance of the error process and allows a better modeling of the properties of the random fluctuations affecting the hemodynamic signal. The practical relevance of our model is illustrated using real event-related fMRI experiments.