code: def beta_covariance(self):
          
          from numpy_sugar.linalg import ddot
  
          tX = self._X["tX"]
          Q = concatenate(self._QS[0], axis=1)
          S0 = self._QS[1]
          D = self.v0 * S0 + self.v1
          D = D.tolist() + [self.v1] * (len(self._y) - len(D))
          D = asarray(D)
          A = inv(tX.T @ (Q @ ddot(1 / D, Q.T @ tX)))
          VT = self._X["VT"]
          H = lstsq(VT, A, rcond=None)[0]
          return lstsq(VT, H.T, rcond=None)[0]
docs: Estimates the covariance-matrix of the optimal beta.
  
          Returns
          -------
          beta-covariance : ndarray
              (Xᵀ(s((1-𝛿)K + 𝛿I))⁻¹X)⁻¹.
  
          References
          ----------
          .. Rencher, A. C., & Schaalje, G. B. (2008). Linear models in statistics. John
             Wiley & Sons.