Existing deep learning-based traffic forecasting models are often trained with MSE as the loss function, which is equivalent to assuming the residual/error to follow independent Gaussian distribution for simplicity of modeling. However, this assumption does not hold especially in traffic forecasting tasks, where the residuals are correlated in both spatial and temporal dimensions. For a multistep forecasting model, we would also expect the variance to increase with the number of steps. In this study, we propose a SpatioTemporal Residual Regularization based on the assumption that the residuals follow a zero-mean multivariate Gaussian distribution with a learnable spatiotemporal covariance matrix. This approach benefits from directly considering correlated spatiotemporal residuals. However, it suffers from scalability issues since the spatiotemporal covariance is often large. For model scalability, we model the spatiotemporal covariance as a sum of Kronecker products of spatial and temporal residual covariance, which significantly reduces the number of parameters and computation complexity. The performance of the proposed method is tested on a traffic speed forecasting task, and the results show that the proposed method improves the model performance by properly dealing with correlated residuals.