A script to compute optimal portfolio weights using both sklearn's Ledoit-Wolf and custom shrinkage estimators

ledoit_wolf_portfolio_rebalancing.py

import yfinance as yf
import numpy as np
from sklearn.covariance import LedoitWolf
from scipy.optimize import minimize
import pandas as pd
from typing import Tuple
from typing import Optional

def shrinkage(returns: np.array) -> Tuple[np.array, float, float]:
    # from https://github.com/WLM1ke/LedoitWolf/
    # Ledoit & Wolf constant correlation unequal variance shrinkage estimator.
    """Shrinks sample covariance matrix towards constant correlation unequal variance matrix.

    Ledoit & Wolf ("Honey, I shrunk the sample covariance matrix", Portfolio Management, 30(2004),
    110-119) optimal asymptotic shrinkage between 0 (sample covariance matrix) and 1 (constant
    sample average correlation unequal sample variance matrix).

    Paper:
    http://www.ledoit.net/honey.pdf

    Matlab code:
    https://www.econ.uzh.ch/dam/jcr:ffffffff-935a-b0d6-ffff-ffffde5e2d4e/covCor.m.zip

    Special thanks to Evgeny Pogrebnyak https://github.com/epogrebnyak

    :param returns:
        t, n - returns of t observations of n shares.
    :return:
        Covariance matrix, sample average correlation, shrinkage.
    """
    t, n = returns.shape
    mean_returns = np.mean(returns, axis=0, keepdims=True)
    returns -= mean_returns
    sample_cov = returns.transpose() @ returns / t

    # sample average correlation
    var = np.diag(sample_cov).reshape(-1, 1)
    sqrt_var = var ** 0.5
    unit_cor_var = sqrt_var * sqrt_var.transpose()
    average_cor = ((sample_cov / unit_cor_var).sum() - n) / n / (n - 1)
    prior = average_cor * unit_cor_var
    np.fill_diagonal(prior, var)

    # pi-hat
    y = returns ** 2
    phi_mat = (y.transpose() @ y) / t - sample_cov ** 2
    phi = phi_mat.sum()

    # rho-hat
    theta_mat = ((returns ** 3).transpose() @ returns) / t - var * sample_cov
    np.fill_diagonal(theta_mat, 0)
    rho = (
        np.diag(phi_mat).sum()
        + average_cor * (1 / sqrt_var @ sqrt_var.transpose() * theta_mat).sum()
    )

    # gamma-hat
    gamma = np.linalg.norm(sample_cov - prior, "fro") ** 2

    # shrinkage constant
    kappa = (phi - rho) / gamma
    shrink = max(0, min(1, kappa / t))

    # estimator
    sigma = shrink * prior + (1 - shrink) * sample_cov

    return sigma, average_cor, shrink


def download_stock_data(stocks: list, start: str, end: str) -> Optional[pd.DataFrame]:
    """
    Download stock data using Yahoo Finance for the provided list of stocks within the given date range.
    """

    print('Initiating stock data download...')
    
    try:
        data = yf.download(stocks, start=start, end=end)
        adj_close_data = data.get('Adj Close')

        if adj_close_data is None:
            print("Adjusted close prices not found in the downloaded data.")
            return None

        print(f"Successfully downloaded data for {adj_close_data.shape[1]} stocks from {start} to {end}.")
        return adj_close_data

    except Exception as e:
        print(f"An error occurred while downloading the stock data: {str(e)}")
        return None


def calculate_returns(data: pd.DataFrame) -> pd.DataFrame:
    """
    Calculate returns for the provided stock data.
    """
    try:
        returns = np.log(data / data.shift(1)).dropna()
        return returns
    except Exception as e:
        print(f"An error occurred while calculating the returns: {e}")
        return None

def calculate_lw_covariance_sklearn(returns: pd.DataFrame) -> pd.DataFrame:
    """
    Calculate the Ledoit-Wolf covariance matrix for the provided returns data.
    Using sklearn's LedoitWolf estimator.
    """
    try:
        lw_cov = LedoitWolf().fit(returns)
        return pd.DataFrame(lw_cov.covariance_, index=returns.columns, columns=returns.columns)
    except Exception as e:
        print(f"An error occurred while calculating the sklearn's Ledoit-Wolf covariance matrix: {e}")
        return None

def calculate_lw_covariance_custom(returns: pd.DataFrame) -> pd.DataFrame:
    """
    Calculate the Ledoit-Wolf covariance matrix for the provided returns data.
    Using the custom shrinkage estimator.
    """
    try:
        sigma, average_cor, shrink = shrinkage(returns.values)
        return pd.DataFrame(sigma, index=returns.columns, columns=returns.columns)
    except Exception as e:
        print(f"An error occurred while calculating the custom Ledoit-Wolf covariance matrix: {e}")
        return None

def portfolio_variance(weights: np.ndarray, cov_matrix: np.ndarray) -> float:
    """
    Define the objective function (portfolio variance)
    """
    return weights @ cov_matrix @ weights


def optimize_portfolio(cov_matrix: np.ndarray, stocks: list) -> np.ndarray:
    """
    Optimize the portfolio to minimize variance.
    """
    # Define the constraints (all weights sum to 1)
    constraints = ({'type': 'eq', 'fun': lambda x: np.sum(x) - 1})
    # Define the bounds for the weights (between 0 and 1)
    bounds = [(0, 1) for _ in range(len(stocks))]
    # Initial guess for the weights
    initial_guess = [1 / len(stocks)] * len(stocks)
    # Optimize!
    try:
        result = minimize(
            portfolio_variance, 
            initial_guess, 
            args=(cov_matrix), 
            bounds=bounds, 
            constraints=constraints)
        return result.x
    except Exception as e:
        print(f"An error occurred while optimizing the portfolio: {e}")
        return None

    
def print_portfolio(stocks: list, weights_sklearn: np.ndarray, weights_custom: np.ndarray):
    """
    Print the optimized portfolio weights for both sklearn and custom covariance estimators.
    """
    print("\n{:<10s} {:<20s} {:<20s}".format('Stock', 'Weight (Sklearn LW)', 'Weight (Custom LW)'))
    print('-' * 50)
    
    for _stock, _weight_sklearn, _weight_custom in zip(stocks, weights_sklearn, weights_custom):
        print("{:<10s} {:<20s} {:<20s}".format(_stock, f'{100*_weight_sklearn:.2f}%', f'{100*_weight_custom:.2f}%'))


def calculate_lw_covariance_sklearn(returns: pd.DataFrame) -> pd.DataFrame:
    """
    Calculate the Ledoit-Wolf covariance matrix for the provided returns data.
    Using sklearn's LedoitWolf estimator.
    """
    try:
        lw_cov = LedoitWolf().fit(returns)
        return pd.DataFrame(lw_cov.covariance_, index=returns.columns, columns=returns.columns)
    except Exception as e:
        print(f"An error occurred while calculating the sklearn's Ledoit-Wolf covariance matrix: {e}")
        return None

    
def calculate_lw_covariance_custom(returns: pd.DataFrame) -> pd.DataFrame:
    """
    Calculate the Ledoit-Wolf covariance matrix for the provided returns data.
    Using the custom shrinkage estimator.
    """
    try:
        sigma, average_cor, shrink = shrinkage(returns.values)
        return pd.DataFrame(sigma, index=returns.columns, columns=returns.columns)
    except Exception as e:
        print(f"An error occurred while calculating the custom Ledoit-Wolf covariance matrix: {e}")
        return None

    
def main():
    # Define the stocks to download.
    stocks = ['AAPL','NFLX','GOOG','AMZN','IBM','INTC','MSFT', 'PACB']
    start = '2022-01-01'
    end = '2023-06-30'

    # Download stock data
    data = download_stock_data(stocks, start, end)

    # Calculate returns
    returns = calculate_returns(data)

    # Compute the Ledoit-Wolf estimator with sklearn
    lw_cov_sklearn = calculate_lw_covariance_sklearn(returns)

    # Optimize the portfolio using sklearn's covariance matrix
    optimal_weights_sklearn = optimize_portfolio(lw_cov_sklearn, stocks)

    # Compute the Ledoit-Wolf estimator with custom shrinkage estimator
    lw_cov_custom = calculate_lw_covariance_custom(returns)

    # Optimize the portfolio using the custom covariance matrix
    optimal_weights_custom = optimize_portfolio(lw_cov_custom, stocks)

    # Print the optimal portfolio weights for both
    print("\nOptimal portfolio weights using sklearn's LedoitWolf and custom shrinkage estimator:")
    print_portfolio(stocks, optimal_weights_sklearn, optimal_weights_custom)


if __name__ == '__main__':
    main()