The authors consider a sequence of \(n\) i.i.d.~observations of a \(p\)-dimensional random vector \((X_i)\), having the factor structure \(X_i=\Lambda\,F_i+\epsilon_i\), where \(\Lambda\) is the loading \(p\times m\) matrix, \(F_i\) is the vector of centred factor variables and \(\epsilon_i\) are the errors -- the idiosyncratic variables. The dimension \(m>0\) is known. The variance is var\((X_i)=\Lambda\,\Lambda'+\Psi\). Factors and idiosyncratic variables are assumed to be uniformly sub-Gaussian. The authors provide \(\ell_1\)-, \(\ell_2\)- and \(\ell_{\infty}\)-error bounds. Their approach is based on a two-step estimation: first the matrices \(\Lambda\) and \(\Psi\) are obtained through a Gaussian quasi-maximum likelihood (QML) estimation (in this step \(\Psi\) is assumed to be diagonal). Conditionally on this first step estimation, the diagonality assumption on \(\Psi\) is relaxed, and by means of various regularisers, both Gaussian QML and least squares loss function are used to obtain a sparse error covariance matrix. The support recovery property is also established. The results are supported by simulations.