Add new shrinkage estimators about pyportfolioopt HOT 7 CLOSED

class LedoitWolfShrinkage:
    """
    Originally described in Ledoit and Wolf (2001).
    The core idea behind Ledoit Wolf Shrinkage is to find a compromise between the sample covariance
    matrix S (although it is the MLE, it has high estimation error) and a shrinkage target F.
    These are then combined with :math:`\delta F + (1-\delta)S`, where :math:`\delta` is
    the shrinkage constant which has to be estimated. The optimal shrinkage constant depends on
    F and S, which is why this module has been designed such that different shrinkage targets are
    sublcasses of LedoitWolfShrinkage.

    https://epublications.bond.edu.au/cgi/viewcontent.cgi?article=1099&context=ejsie

    The default is constant correlation shrinkage, as per Ledoit and Wolf 2004.
    """

    # TODO implement a default identity shrinkage

    def __init__(self, prices):
        self.X = prices.pct_change().dropna(how="all")
        self.S = self.X.cov()  # sample cov is the MLE
        self.F = None  # default shrinkage target
        self.delta = None  # shrinkage constant

    def _shrinkage_target(self):
        """
        Calculates the constant correlation shrinkage target, as per Ledoit and Wolf 2004.
        All pairwise correlations are assumed to be equal to the mean of all correlations,
        and the diagonal is just the std of each asset's returns.
        Note that these functions particularly depend on X being an pandas dataframe, because
        it provides functions that are very robust to missing data.
        :return: constant correlation shrinkage target
        :rtype: np.array
        """
        T, N = self.X.shape
        correlation_matrix = self.X.corr()
        # r_bar is the average pairwise correlation (excluding diagonal)
        r_bar = (np.sum(correlation_matrix.values) - N) / (N * (N - 1))
        # r_bar = (np.mean(correlation_matrix.values)  average correlation

        # Use numpy's broadcasting to calculate s_ii * s_jj
        # Then f_ij = r_bar * sqrt(s_ii * s_jj)
        F = r_bar * np.sqrt(self.X.std().values *
                            self.X.std().values.reshape((-1, 1)))
        # f_ii = s_ii
        np.fill_diagonal(F, self.X.std().values)
        self.F = F
        return self.F, r_bar

    def _pi_matrix(self):
        """
        Estimates the asymptotic variances of the scaled entries of the sample covariance matrix.
        This is used to calculate pi and rho.
        :return: scaled asymptotic variances, in an NxN array where N is the number of assets.
        :rtype: np.array
        """

        T, N = self.X.shape
        Xc = np.array(self.X - self.X.mean())  # centre X
        Xc = np.nan_to_num(Xc)
        S = np.array(self.S)
        pi_matrix = np.zeros((N, N))
        for i in range(N):
            for j in range(i, N):
                pi_matrix[i, j] = pi_matrix[j, i] = np.sum(
                    (Xc[:, i] * Xc[:, j] - S[i, j]) ** 2
                )
        return pi_matrix / T

    def _pi(self):
        """
        Estimates the sum of asymptotic variances of the scaled entries of the
        sample covariance matrix. One of the three components of the optimal shrinkage
        constant
        :return: pi(hat)
        :rtype: float
        """
        pi_matrix = self._pi_matrix()
        return np.sum(pi_matrix)

    # TODO implement
    def _rho(self):
        Xc = self.X - self.X.mean()
        T, N = self.X.shape

        pi_diag = np.trace(self._pi_matrix())
        _, r_bar = self._shrinkage_target()

        asycov = 0

        for i in range(N):
            for j in range(N):
                if j == i:
                    continue

                theta_ii = np.sum(
                    (Xc.iloc[:, i] ** 2 - self.S.iloc[i, i])
                    * (Xc.iloc[:, i] * Xc.iloc[:, j] - self.S.iloc[i, j])
                )
                theta_jj = np.sum(
                    (Xc.iloc[:, j] ** 2 - self.S.iloc[j, j])
                    * (Xc.iloc[:, i] * Xc.iloc[:, j] - self.S.iloc[i, j])
                )

                first_term = np.sqrt(
                    self.S.iloc[j, j] / self.S.iloc[i, i]) * theta_ii
                second_term = np.sqrt(
                    self.S.iloc[i, i] / self.S.iloc[j, j]) * theta_jj
                asycov += first_term + second_term

        return pi_diag + r_bar / 2 * asycov / T

    def _gamma(self):
        """
        Gamma estimates deviance of the shrinkage target (i.e sum of element-wise squared distances
        between shrinkage target matrix and sample covariance matrix).
        :return: gamma
        :rtype: float
        """
        if self.F is None:
            self._shrinkage_target()
        S = np.array(self.S)
        return np.sum((self.F - S) ** 2)

    def _optimal_shrinkage_constant(self):
        """
        Calculates the optimal value of the shrinkage constant, as per Ledoit and Wolf. This
        is truncated to [0, 1], in the rare case that k/T is outside of this range.
        :return: optimal shrinkage constant
        :rtype: float between 0 and 1.
        """
        T = self.X.shape[0]
        pi = self._pi()
        rho = self._rho()
        gamma = self._gamma()
        kappa = (pi - rho) / gamma
        self.delta = max(0, min(kappa / T, 1))
        return self.delta

    def shrink(self):
        if self.delta is None:
            self._optimal_shrinkage_constant()
        # Return shrunk matrix
        S_hat = self.delta * self.F + (1 - self.delta) * self.S
        return S_hat