Accurate covariance matrices are required for a reliable estimation of cosmological parameters from pseudo-power spectrum estimators. In this work, we focus on the analytical calculation of covariance matrices. We consider the case of observations of the Cosmic Microwave Background in temperature and polarization on a small footprint such as in the SPT-3G experiment, which observes 4% of the sky. Power spectra evaluated on small footprints are expected to have large correlations between modes, and these need to be accurately modelled. We present, for the first time, an algorithm that allows an efficient (but computationally expensive) exact calculation of analytic covariance matrices. Using it as our reference, we test the accuracy of existing fast approximations of the covariance matrix. We find that, when the power spectrum is binned in wide bandpowers, current approaches are correct up to the 5% level on the SPT-3G small sky footprint. Furthermore, we propose a new approximation which improves over the previous ones reaching a precision of 1% in the wide bandpowers case and generally more than 4 times more accurate than current approaches. Finally, we derive the covariance matrices for mask-corrected power spectra estimated by the PolSpice code. In particular, we include, in the case of a small sky fraction, the effect of the apodization of the large scale modes. While we considered the specific case of the CMB, our results are applicable to any other cosmological probe which requires the calculation of pseudo-power spectrum covariance matrices.