MVO and EW MVO

We’ve spilled enough pixelated ink to forego a lengthy discussion of the MVO portfolio algorithm. But suffice it to say, using convex optimization, one finds the weights that generate the maximum return for a given level of risk. The formula is:

\(\text{max: }\mu^T w – \gamma w^T\Sigma w\)
\(\text{subject to: } 1^Tw = 1\)

We won’t go into too much detail on the equation, as we’ve dealt with it before. While this is a general formula with \(\gamma\) representing a risk aversion metric, for this post we’ll actually be maximizing the Sharpe ratio or:

\(\text{max: }\frac{\mu^Tw}{{(w^T\Sigma w)}^{1/2}}\)

We’ll also use an exponentially-weighed version of returns to estimate future returns for the MVO.

Kelly criterion

Not usually used to find portfolio weights, the Kelly Criterion is discussed frequently by fundamental investors. Developed by John L. Kelly Jr. at Bell Labs to parse noise from signal, it found success at Vegas blackjack tables until the bosses figured it out, broke some knuckles, and added more decks to the dealer’s hand. The idea is to ask, how much should one “bet” as percent of his or her bankroll given an expected edge and desire never to go bankrupt. The elegantly simple formula is:

\(K\% = w – \frac{1-w}{r}\)

Where \(w\) is your expected win rate and \(r\) is the ratio of your win rate to loss rate. Note that if you’re expected win rate is 50%, the formula returns a 0 weight, or no bet. Smart formula.

We create an asset level Kelly, by deriving the win rate as the frequency of positive returns for the series. We calculate that for every asset and then scale the asset Kelly by the sum of all the asset Kellys to get a portfolio weight for each asset like the following:

\(\text{Kelly weight: }\frac{k_{i}}{\sum_{i=1}^{n} k_{i}}\)

We could get a little more sophisticated, but it’s been a long week.

Black-Litterman Model

Created by Fischer Black of Black-Scholes option pricing fame and Robert Litterman, the model seeks to improve on MVO. While it is impossible to do the model justice, think of it this way. The problem with MVO is that if you get your expected return estimate wrong you’re screwed. Mess up the estimates on your covariance matrix? That’s bad too, and there are a bunch of techniques to handle that. But bad return estimates tend to dominate the results². Enter Black-Litterman, which tries to solve the return estimation problem by applying Bayesian methods³ to update one’s views on returns given the market price of risk.

Here’s the logic. Investors have views on returns. Probably not overly accurate views. But they’re more likely to be right about whether a set of assets will outperform another and which ones might enjoy higher excess returns. The model incorporates these views, by updating the market implied excess return base rate (the prior) with the forecasts modified by the confidence of that forecast, to return a new estimate of expected excess returns (the posterior). Here’s the formula.

\(E[R] = [(\tau\Sigma)^{-1} + P^T\Omega^{-1}P]^{-1}[(\tau\Sigma)^{-1}\Pi + P^T\Omega^{-1}Q]\)

Don’t ask us to explain it! Essentially, it says expected returns (\(E[R]\)) are a function of prior market-implied expected returns (\(\Pi\)) and the investor’s views (\(Q\)) adjusted for uncertainty (\(\Omega\)) on how much to weight (\(P\)) a particular asset. The covariance of the assets is also in there (\(\Sigma\)) along with another weighting mechanism (\(\tau\)) that dictates the amount by which one should emphasize active vs. passive risk.⁴ Clear?

The important point about the BL model is that it’s not an optimization algorithm. It’s a way to get better return estimates to put into an MVO or other algorithm.

Hierarchical Risk Parity

This model is the brainchild of Marcos López de Prado, Cornell professor, Quant of the Year, and author of Advances in Financial Machine Learning. Again, very hard to encapsulate in a short paragraph. The main idea is to solve the portfolio allocation decision using machine learning instead of classical analytical techniques. At a basic level. Prado argues that MVO is unreliable because its leads to unstable solutions, asset concentrations that expose the portfolio to idiosyncratic risk, and out-of-sample underperformance. The reasons this occurs has to do with the math.⁵ But the intuition is when weights are unconstrained and without context, that can lead to substantially different solutions even when making small changes to the estimates. Prado solves this by introducing a hierarchical structure that prevents too much exposure to idiosyncratic risk.

Think of it this way. The traditional MVO solution would treat Exxon, Chevron, and ConocoPhillips as independent and might very well assign high weights to them in a high and rising oil environment. But a human PM would be worried about having too much energy exposure, especially when oil (inevitably) tanks. HRP is supposed to solve this, producing better diversification out-of-sample.⁶

Examining the models

So what do all these portfolio allocation algorithms look like? First, let’s visualize our starting point: the cumulative return graph of our four asset classes during the base period to run the algorithms: 1987-1991.

Stocks had a great run even if the ride was bit bumpy. Bonds weren’t bad either (Oh the days when bonds actually yielded something!). Real estate ticked along nicely, while the gold was a disappointment. Given these results, what weights would the various algorithms have recommended? Recall, for the MVO and BL portfolios, we’ve set the algorithm to pick the portfolio with the highest Sharpe ratio. The BL model requires us to input our view on expected returns and our confidence on those views, rather than interrupt the narrative, the reader can find these inputs after the code.

Interestingly, most of the portfolios select a high weighting to real estate. I don’t think any fundamental investor would have chosen the same weight. The Kelly portfolio seems to be the most reasonable. The exponentially-weighted MVO (EW MVO) features the lowest weighting to real estate after the naive portfolio. Although the weights on the naive portfolio should be obvious, we include them as reference.

Now let’s look at how these portfolios perform over the next period: 1992-1996. First, here’s the cumulative return graph.

It should be relatively obvious that any portfolio that overweights the stock market relative to others will perform better. We show two versions of the weightings below. One with no rebalancing, the other with a yearly rebalance back to the original scheme. Note, that no rebalancing allows the weights to float, leaving the ending exposure often very different to the initial one.

As is apparent, the Naive and Kelly portfolios outperform the more sophisticated algorithms given their higher exposure to stocks; though EW MVO wasn’t bad. We don’t think we need to run a t-test to prove that this outperformance was significant at around 400bps. One more thing to note is that yearly rebalancing underperformed no rebalancing. Something to consider especially if your robo-advisor advocates frequent rebalancing. Whatever the case, how did the portfolios perform in terms of the Sharpe ratio?

Not as well.⁷ The sophisticated algorithms clearly outperform the less sophisticated ones, except for EW MVO, which surprisingly performs worse than the Kelly portfolio. This sort of proves why you don’t want to rely too much on recent performance for your asset weightings.

Given the return outperformance and Sharpe underperformance of the simple weighting schemes relative to the sophisticated ones, it should be obvious that volatility is the main reason for the difference. We show the “vol” comparisons below.

As they should, the algorithms meant to minimize risk for a given level of return, do exactly what they say on the tin. The simple algorithms clearly have significantly higher volatility on relative basis, but the overall base is quite low.

The key questions begged by these results are the following. Is it rational to suffer the slings and arrows of lower Sharpe ratios to enjoy a huge return advantage? Normative finance says no. But great investors say yes. Warren Buffett famously said he’d be happy with a volatile 15% return over a smooth 12% one.

Second, why did the Black-Litterman model and the Hierarchical Risk Parity perform similarly? They’re different approaches with different algorithms. This might be an artifact of the data or the expected return and conviction we used for the BL model. We suspect part of the reason is the paucity of assets. If there were more assets, the hierarchical algorithm would probably have differentiated itself as there would be, well, more hierarchies to distinguish between.

Third, is the Naive ratio, so naive after all? If we don’t really know the future, why not keep our asset weights on an even keel so we never miss the next wave? Why not minimize the probability of future regret as one Nobel Prize winner did?

We’ll end it here given that we’ve consumed more pixels than planned. In our next post, we’ll iterate through the other periods to get a total period performance of the different weighting schemes. If you have questions, you can reach us at the email address below. Until next time, here’s the code:

Built using R 4.0.3, and Python 3.8.8

[R]
## Load R libraries
suppressPackageStartupMessages({
  library(tidyquant)
  library(tidyverse)
  library(reticulate)
})

# Apologies to R users. This post is all Python. If anyone knows of R library that approximates PyPortfolioOpt, please let us know!

[Python]
## Load packages
import pandas as pd
import pandas_datareader.data as dr
import numpy as np
import yfinance as yf
import pandas_datareader.famafrench as ffr
import cvxpy as cp
import matplotlib.pyplot as plt

plt.rcParams['figure.figsize'] = (12,6)
plt.style.use('ggplot')

## Save graphs for blog
import os
DIR = "C:/Users/user/path/to/blog/static" 

def save_fig_blog(fig_id, tight_layout=True, fig_extension="png", resolution=300):
    path = os.path.join(DIR, fig_id + "." + fig_extension)
    print("Saving figure", fig_id)
    if tight_layout:
        plt.tight_layout()
    plt.savefig(path, format=fig_extension, dpi=resolution)
    
## Create dataset
start_date = '1970-01-01'
end_date = '2020-12-31'
symbols = ["WILL5000INDFC", "BAMLCC0A0CMTRIV", "GOLDPMGBD228NLBM", "CSUSHPINSA"]
sym_names = ["stock", "bond", "gold", "realt"]
filename = 'port_const_25.pkl'

try:
    df = pd.read_pickle(filename)
    print('Data loaded')
except FileNotFoundError:
    print("File not found")
    print("Loading data", 30*"-")
    data = dr.DataReader(symbols, 'fred', start_date, end_date)
    data.columns = sym_names
    data_mon = data.resample('M').last()
    
    dataf = data_mon.pct_change()
    
    # Pull Fama-French risk-free rate. This series is different than what we find on FRED. We figure FF know what they're doing, so use their data!
    ff_d = ffr.FamaFrenchReader('F-F_Research_Data_Factors_daily', start=start_date, end=end_date)
    ff = ff_d.read()[0]
    ff_mon = ff.resample('M').last()
    df = pd.merge(dataf, ff_mon, how= 'left', left_index=True, right_index=True)
    df.drop(['Mkt-RF','SMB','HML'], axis=1, inplace=True)
    df.rename(columns= {'RF' : 'rfr'}, inplace=True)
    df['rfr'] = new_df['rfr']/100 

# Save for later use
df.loc['1987':].to_pickle('port_const_25.pkl')

## Create first series
df_1 = df.loc['1987':'1991', 'stock':'realt']

# Plot cumulative returns
(((1+df_1).cumprod()-1)*100).plot(cmap="Blues_r")
plt.xlabel("")
plt.ylabel("Return (%)")
plt.legend(['Stocks', 'Bonds', 'Gold', 'Real estate'], loc='upper center', ncol=4)
plt.title('Cumulative returns to the four asset classes: 1987-1991')
plt.show()

## MVO
# Function to return weights for maximum Sharpe ratio (tangency) portfolio
def max_sharpe(mu, sigma, samples=100):
    w = cp.Variable(len(mu))
    gamma = cp.Parameter(nonneg=True)
    ret = mu.T @ w 
    risk = cp.quad_form(w, sigma)
    objective = cp.Maximize(ret - gamma*risk)
    prob = cp.Problem(objective, [cp.sum(w) == 1, w >= 0])
    
    risk_data = np.zeros(samples)
    ret_data = np.zeros(samples)
    wt_data = np.zeros(samples, dtype=object)
    gamma_vals = np.logspace(-2, 3, num=samples)
    for i in range(samples):
        gamma.value = gamma_vals[i]
        prob.solve()
        risk_data[i] = cp.sqrt(risk).value
        ret_data[i] = ret.value
        wt_data[i] = w.value
    
    return wt_data[(ret_data/risk_data).argmax()]

# Run function    
mu = df_1.mean().values # CVXPY needs numpy arrays
cov = df_1.cov().values
eff_sharpe = max_sharpe(mu, cov)

## EW MVO
mu_ema = df_1.ewm(12).mean().iloc[-1].values
ema_sharpe = max_sharpe(mu_ema, cov)

## Kelly Criterion
# Create function
def kelly(series):
    win = sum(series > 0)/len(series)
    loss = 1 - win
    kelly = win - loss/(win/(loss+1e-6)) 
    return  0 if kelly < 0 else kelly 

# Run function    
kelly_wts = df_1.apply(kelly)/df_1.apply(kelly).sum()

## Black-Litterman
from pypfopt import black_litterman, risk_models
from pypfopt import BlackLittermanModel, plotting

# Create function
# We're sorry we didn't go into more detail on how this works. 
def get_bl_wts(return_data, mkt_cap_wt = np.array([0.12,0.48,0.5,0.35]), 
               views = {'stock':0.07, 'bond':0.05, 'gold':0.0, 'realt':0.04},
               confidence = [0.75, 0.65, 0.5, 0.55]):
    delta = (return_data.mean())/(return_data.var())
    cov_mat = return_data.cov()
    mkt_prior = delta * (cov_mat.values @ mkt_cap_wt)

    bl = BlackLittermanModel(cov_mat, pi = mkt_prior,
                        absolute_views = views,
                        omega='idzorek',
                        view_confidences=confidence)
    
    ret_bl = bl.bl_returns()
    cov_bl = bl.bl_cov()
    return max_sharpe(ret_bl.values, cov_bl.values)

## Hierarchical Risk Parity
from pypfopt import HRPOpt

def get_hrp_wts(return_data):
    hrp = HRPOpt(return_data)
    hrp.optimize()
    return pd.Series(hrp.clean_weights())

## Graph weight comparisons
def wt_comp_plot(array, col_names = ['MVO', 'EW MVO', 'Kelly', 'BL'], ind_names = ['Stock', 'Bond', 'Gold', 'Real Estate'],
                 save_fig = False, png_name = None):
    (pd.DataFrame(array*100, 
            columns = col_names,
            index = ind_names)
            .plot(kind='bar', cmap="Blues_r", rot=0))
    plt.legend(loc='upper center', ncol=len(col_names))
    plt.ylabel("Weight (%)")
    plt.title("Asset weights by optimization method")
    if save_fig:
        save_fig_blog(png_name)
    plt.show()

## Test the optimizations

# create test set    
for_weights = np.c_[weights]
df_2 = df.loc['1992':'1996', 'stock':'realt']
rfr_2 = df.loc['1992':'1996', 'rfr']

# Create rebalancing function
def rebal_func(act_ret, weights, frequency=None):
    """Rebalancing logic: 
       If not rebalancing, as value of assets change, original weight in portfolio changes.
       Need to adjust the weight for the overall return in the portfolio to calculate next period return.
       Hence (weight_asset_i + return_asset_i) / (1 + return_portfolio).
       If rebalancing, just multiply be original weight."""
    
    beg_weights = weights
    ret_vec = np.zeros(len(act_ret))
    wt_mat = np.zeros((len(act_ret), len(act_ret.columns)))
    for i in range(len(act_ret)):
        if frequency:
            if (i % frequency == 0) | (i < 1):
                wt_ret = act_ret.iloc[i,:].values*beg_weights
            else:
                wt_ret = act_ret.iloc[i,:].values*new_weights
        else:
            if i < 1 :
                wt_ret = act_ret.iloc[i,:].values*beg_weights
            else:
                wt_ret = act_ret.iloc[i,:].values*new_weights
            
        ret = np.sum(wt_ret)
        ret_vec[i] = ret
        
        if frequency:
            if (i % frequency == 0) | (i < 1):
                new_weights = (beg_weights + wt_ret)/(1 + ret)

            else:
                new_weights = (new_weights + wt_ret)/(1 + ret)
        else:
            if i < 1 :
                new_weights = (beg_weights + wt_ret)/(1 + ret)
            else:
                new_weights = (new_weights + wt_ret)/(1 + ret)
            
        wt_mat[i,] = new_weights
    
    return ret_vec, wt_mat
    
# Create function to port rebalancing results into data frame    
def get_rebal_df(dataf, port_weights, ports=['Naive','MVO', 'EW MVO', 'Kelly', 'BL', 'HRP'], frequency = None):
    if (port_weights.shape[1] != len(ports)) | (port_weights.shape[0] != len(dataf.columns)):
        print("Weights and portfolios don't match!")
    else:
        results = {}

        for i in range(len(ports)):
            rets, wts = rebal_func(dataf, port_weights.T[i], frequency=frequency)
            results[ports[i]] = {'returns':rets, 'weights': wts}
            
        port_rets = []
        for key in results.keys():
            rets = results[key]['returns'].mean()*12
            risk = results[key]['returns'].std()*np.sqrt(12)
            sharpe = (rets)/risk
            port_rets.append([rets, risk, sharpe])

        return pd.DataFrame(port_rets, columns=['return', 'risk', 'sharpe'], index = ports)

port_df = get_rebal_df(df_2, for_weights)
port_rebal_df = get_rebal_df(df_2, for_weights, frequency=12)

# Create function plot results
def port_comp_plot(dataf1, metric, dataf2 = None, ports = ['Naive','MVO', 'EW MVO', 'Kelly', 'BL', 'HRP'],
                   save_fig=False, png_name=None):
    if dataf2 is not None:
        (pd.DataFrame(np.c_[dataf1[metric], dataf2[metric]], 
                 columns = ['None', 'Yearly'], 
                 index = ports)
                 .plot(kind='bar', color=['blue', 'grey'], rot=0))
    else:
        (pd.DataFrame(np.c_[dataf1[metric]], 
                 index = ports)
                 .plot(kind='bar', color='blue', rot=0, legend=None))
    
    metric_dict = {'return':'Return', 'sharpe':'Sharpe ratio', 'risk': 'Volatility'}
    metric_label = {'return':'Percent (%)', 'sharpe':'Ratio', 'risk': 'Percent (%)'}
    
    plt.ylabel(metric_label[metric])
    
    if dataf2 is not None:
        plt.title("{} comparisons by portfolio and rebalancing regime".format(metric_dict[metric]))
    else:
        plt.title("{} comparisons by portfolio".format(metric_dict[metric]))
        
    if save_fig:
        save_fig_blog(png_name)
    
    plt.show()
    
# Plot results: returns, Sharpe ratio, volatility    
port_comp_plot(port_df*100, 'return', port_rebal_df*100)

port_comp_plot(port_df, 'sharpe', port_rebal_df)

port_comp_plot(port_df*100, 'risk', port_rebal_df*100)