Accurately modeling the correlation structure of errors is essential for reliable uncertainty quantification in probabilistic time series forecasting. Recent deep learning models for multivariate time series have developed efficient parameterizations for time-varying contemporaneous covariance, but they often assume temporal independence of errors for simplicity. However, real-world data frequently exhibit significant error autocorrelation and cross-lag correlation due to factors such as missing covariates. In this paper, we present a plug-and-play method that learns the covariance structure of errors over multiple steps for autoregressive models with Gaussian-distributed errors. To achieve scalable inference and computational efficiency, we model the contemporaneous covariance using a low-rank-plus-diagonal parameterization and characterize cross-covariance through a group of independent latent temporal processes. The learned covariance matrix can be used to calibrate predictions based on observed residuals. We evaluate our method on probabilistic models built on RNN and Transformer architectures, and the results confirm the effectiveness of our approach in enhancing predictive accuracy and uncertainty quantification without significantly increasing the parameter size.