def cov_random(self, D, Sinv=None):
        if Sinv is not None:
            self.compute_P(Sinv)
        t = np.dot(self.Z, D)
        return D - np.dot(np.dot(t.T, self.P), t)
Approximate covariance of estimates of random effects. Just after
dataset/ETHPy150Open statsmodels/statsmodels/statsmodels/sandbox/panel/mixed.py/Unit.cov_random