Calculate the Constant-Correlation covariance matrix.

        :return: shrunk sample covariance matrix
        :rtype: np.ndarray

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)