I am trying to create a Simulated Portfolio Optimization based on Efficient Frontier on 50 stocks, which you can find the csv here. Yet it already takes me several minutes to get a suboptimal solution: I can't draw an accurate efficient frontier:

enter image description here

Whereas it should be something like:

enter image description here

So basically, I want to create an efficient frontier of optimization of the weights $w_i$ of stocks in a portfolio of actions $i$ which returns are $x_i$.

I have imagined there is another way to get the weights in the following way. It should be easier to get that efficient frontier getting weights with given, fixed, portfolio standard deviations $\sigma_p$. Indeed, one can fix a grid of volatilities $σ_{p_1},...σ_{p_n}$, then for each $σ_{p_i}$, maximize expected returns with the constraint that the volatility is no larger than $σ_{p_i}$, to get $μ_{p_i}$. Then $(σ_{p_i},μ_{p_i})$ are $n$ points on the efficient frontier.

So, a first step would be to get the weights for one volatility $σ_{p}$. Knowing that for two assets, the portfolio variance $\sigma_p$ is

$$ \begin{align} \sigma_p &= \sqrt{w_1^2\sigma_1^2 + w_2^2\sigma_2^2 +2w_1w_2cov(x_1,x_2)}\\ \end{align} $$

Where $\forall i\neq p,\sigma_i$ are the standard deviations for a given asset.

We can maximize returns $r$ which are equals to the weights time the individual results for each action $RW$. This leads to the following optimization problem (I reduced it to two variables for the sake of simplicity):

$$\begin{cases}\max r\\ &\sigma_p \leq value\\ &\sigma_p = \sqrt{w_1^2\sigma_1^2+w_2^2\sigma^2+2w_1w_2cov_{1,2}}\\ &r = w_1r_1+w_2r_2\\ &\forall i, w_i\geq 0 \end{cases}$$

I don't know how to write it in matrix formulation:

$$\begin{cases}\max r\\ &\sigma_p \leq value\\ &\sigma_p = \sqrt{W^2\Sigma^2+2WW^TCOV}\\ &r = WR\\ &\forall i, w_i\geq 0 \end{cases}$$

Where COV is the covariance matrix between all assets.

But I don't know if it's right and how to write it in python.


My original approach was naive sampling. It doesn't work well because the efficient frontier is a very small subspace of the space I'm exploring:

import pandas as pd  
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import quandl
#import scipy.optimize as scoplt.style.use('fivethirtyeight')

def portfolio_annualised_performance(weights, mean_returns, cov_matrix):
    returns = np.sum(mean_returns*weights ) *252
    std = np.sqrt(np.dot(weights.T, np.dot(cov_matrix, weights))) * np.sqrt(252)
    return std, returns

def random_portfolios(num_portfolios, mean_returns, cov_matrix, risk_free_rate, df):
    results = np.zeros((3,num_portfolios))
    weights_record = []
    for i in range(num_portfolios):
        weights = np.random.random(len(df.columns))
        weights /= np.sum(weights)
        portfolio_std_dev, portfolio_return = portfolio_annualised_performance(weights, mean_returns, cov_matrix)
        results[0,i] = portfolio_std_dev
        results[1,i] = portfolio_return
        results[2,i] = (portfolio_return - risk_free_rate) / portfolio_std_dev
    return results, weights_record

def display_simulated_ef_with_random(mean_returns, cov_matrix, num_portfolios, risk_free_rate, df):
    results, weights = random_portfolios(num_portfolios,mean_returns, cov_matrix, risk_free_rate, df)

    max_sharpe_idx = np.argmax(results[2])
    sdp, rp = results[0,max_sharpe_idx], results[1,max_sharpe_idx]
    print("results[0,max_sharpe_idx], results[1,max_sharpe_idx]: ", results[0,max_sharpe_idx], results[1,max_sharpe_idx])
    max_sharpe_allocation = pd.DataFrame(weights[max_sharpe_idx],index=df.columns,columns=['allocation'])
    max_sharpe_allocation.allocation = [round(i*100,2)for i in max_sharpe_allocation.allocation]
    max_sharpe_allocation = max_sharpe_allocation.T

    min_vol_idx = np.argmin(results[0])
    sdp_min, rp_min = results[0,min_vol_idx], results[1,min_vol_idx]
    min_vol_allocation = pd.DataFrame(weights[min_vol_idx],index=df.columns,columns=['allocation'])
    min_vol_allocation.allocation = [round(i*100,2)for i in min_vol_allocation.allocation]
    min_vol_allocation = min_vol_allocation.T

    print("Maximum Sharpe Ratio Portfolio Allocation\n")
    print("Annualised Return:", round(rp,2))
    print("Annualised Volatility:", round(sdp,2))
    print("Minimum Volatility Portfolio Allocation\n")
    print("Annualised Return:", round(rp_min,2))
    print("Annualised Volatility:", round(sdp_min,2))

    plt.figure(figsize=(10, 7))
    plt.scatter(results[0,:],results[1,:],c=results[2,:],cmap='YlGnBu', marker='o', s=10, alpha=0.3)
    plt.scatter(sdp,rp,marker='*',color='r',s=500, label='Maximum Sharpe ratio')
    plt.scatter(sdp_min,rp_min,marker='*',color='g',s=500, label='Minimum volatility')
    plt.title('Simulated Portfolio Optimization based on Efficient Frontier')
    plt.xlabel('annualised volatility')
    plt.ylabel('annualised returns')

    return max_sharpe_allocation, min_vol_allocation

returns = df.pct_change()
mean_returns = returns.mean()
cov_matrix = returns.cov()
num_portfolios = 750000
risk_free_rate = 0.0178

min_vol_al, max_sharpe_al = display_simulated_ef_with_random(mean_returns, cov_matrix, num_portfolios, risk_free_rate, df)

As a side note, one should also notice that:

std = np.sqrt(np.dot(weights.T, np.dot(cov_matrix, weights))) * np.sqrt(252)

Which leads to another equation which I don't know if it can be useful:

$$W^TCW = (\frac{\sigma_p}{\sqrt{252}})^2$$

So I wonder if we can rewrite:

$$\begin{align} \sigma_p &= \sqrt{W^2\Sigma^2+2WW^TCOV}\\ \Leftrightarrow \sigma_p &= \sqrt{W^2\Sigma^2+2(\frac{\sigma_p}{\sqrt{252}})^2}\\ \Leftrightarrow W^2\Sigma^2 &= \sigma_p^2 -2(\frac{\sigma_p}{\sqrt{252}})^2\\ \Leftrightarrow W &= \sqrt{\frac{\sigma_p^2 -2(\frac{\sigma_p}{\sqrt{252}})^2}{\Sigma^2}} \end{align}$$

But when trying this out with Python and with $\sigma_p=0.2$:

W = np.sqrt((0.2**2-2*(0.2/np.sqrt(252)**2))/cov_matrix)

It rather looks like a covariance matrix:

enter image description here

  • $\begingroup$ is this for fun or work? $\endgroup$ Jun 3, 2020 at 13:47
  • $\begingroup$ @Aksakalalmostsurelybinary Hmm, let's say both! It's for personal investments, and learning as well. $\endgroup$ Jun 3, 2020 at 13:51
  • $\begingroup$ if it's for work, then read the discussion of Eq(30) here the point is that the ordinary estimate of covaraince matrix "as is" is not a robust tool to construct a portfolio. that's one reason why we have approaches like Black-Litterman $\endgroup$ Jun 3, 2020 at 13:59
  • $\begingroup$ @Aksakalalmostsurelybinary Thanks! Wow, that seems too inteligent for me. So I'm going to say it's rather for fun! Interesting, I'm living nearby the address they mention at the top of the paper :p $\endgroup$ Jun 3, 2020 at 14:02
  • $\begingroup$ it's actually not as scary as it may seem at first sight. consider this: you obtained some daily returns, then calculated the covariance matrix. your data set was random, i.e. if you obtain the same series in different period you get a different set of values, and (!) a different covariance matrix. hence, your covariance matrix is itself a random object. you're lucky if it's a good estimate of a true covariance matrix. now, when you build your portfolio from this matrix as in textbooks, they rarely tell you that the results is most likely a garbage if the set of stock is large. $\endgroup$ Jun 3, 2020 at 14:06


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.