We consider the estimation of approximate factor models for time series data, where strong serial and cross-sectional correlations amongst the idiosyncratic component are present. This setting comes up naturally in many applications, but existing approaches in the literature rely on the assumption that such correlations are weak, leading to mis-specification of the number of factors selected and consequently inaccurate inference. In this paper, we explicitly incorporate the dependent structure present in the idiosyncratic component through lagged values of the observed multivariate time series. We formulate a constrained optimization problem to estimate the factor space and the transition matrices of the lagged values {\em simultaneously}, wherein the constraints reflect the low rank nature of the common factors and the sparsity of the transition matrices. We establish theoretical properties of the obtained estimates, and introduce an easy-to-implement computational procedure for empirical work. The performance of the model and the implementation procedure is evaluated on synthetic data and compared with competing approaches, and further illustrated on a data set involving weekly log-returns of 75 US large financial institutions for the 2001-2016 period.