docs: def shrink(self):
          """
          Calculate the Constant-Correlation covariance matrix.
  
          :return: shrunk sample covariance matrix
          :rtype: np.ndarray
          """
          x = np.nan_to_num(self.X.values)
  
          # de-mean returns
          t, n = np.shape(x)
          meanx = x.mean(axis=0)
          x = x - np.tile(meanx, (t, 1))
          xmkt = x.mean(axis=1).reshape(t, 1)
  
          # compute sample covariance matrix
          sample = np.cov(np.append(x, xmkt, axis=1), rowvar=False) * (t - 1) / t
          covmkt = sample[0:n, n].reshape(n, 1)
          varmkt = sample[n, n]
          sample = sample[:n, :n]
          prior = np.dot(covmkt, covmkt.T) / varmkt
          prior[np.eye(n) == 1] = np.diag(sample)
  
          # compute shrinkage parameters
          if self.delta is None:
              c = np.linalg.norm(sample - prior, "fro") ** 2
              y = x ** 2
              p = 1 / t * np.sum(np.dot(y.T, y)) - np.sum(sample ** 2)
              # r is divided into diagonal
              # and off-diagonal terms, and the off-diagonal term
              # is itself divided into smaller terms
              rdiag = 1 / t * np.sum(y ** 2) - sum(np.diag(sample) ** 2)
              z = x * np.tile(xmkt, (n,))
              v1 = 1 / t * np.dot(y.T, z) - np.tile(covmkt, (n,)) * sample
              roff1 = (
                  np.sum(v1 * np.tile(covmkt, (n,)).T) / varmkt
                  - np.sum(np.diag(v1) * covmkt.T) / varmkt
              )
              v3 = 1 / t * np.dot(z.T, z) - varmkt * sample
              roff3 = (
                  np.sum(v3 * np.dot(covmkt, covmkt.T)) / varmkt ** 2
                  - np.sum(np.diag(v3).reshape(-1, 1) * covmkt ** 2) / varmkt ** 2
              )
              roff = 2 * roff1 - roff3
              r = rdiag + roff
  
              # compute shrinkage constant
              k = (p - r) / c
              shrinkage = max(0, min(1, k / t))
              self.delta = shrinkage
          else:
              # use specified constant
              shrinkage = self.delta
  
          # compute the estimator
          sigma = shrinkage * prior + (1 - shrinkage) * sample
          return self.format_and_annualise(sigma)