1
$\begingroup$

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.

Context

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')
np.random.seed(777) 

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)
        weights_record.append(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("-"*80)
    print("Maximum Sharpe Ratio Portfolio Allocation\n")
    print("Annualised Return:", round(rp,2))
    print("Annualised Volatility:", round(sdp,2))
    print("\n")
    print(max_sharpe_allocation)
    print("-"*80)
    print("Minimum Volatility Portfolio Allocation\n")
    print("Annualised Return:", round(rp_min,2))
    print("Annualised Volatility:", round(sdp_min,2))
    print("\n")
    print(min_vol_allocation)

    plt.figure(figsize=(10, 7))
    plt.scatter(results[0,:],results[1,:],c=results[2,:],cmap='YlGnBu', marker='o', s=10, alpha=0.3)
    plt.colorbar()
    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')
    plt.legend(labelspacing=0.8)

    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

Attempt with LinAlg's answer

I think LinAlg's answer is great as it uses the fact that the minimizer of a quadratic objective under linear constraints can be obtained by solving a linear system. Which is really easy to compute and not that hard to understand. To get the best sharpe ratio it gives the following:

def efficient_portfolios(returns, risk_free_rate, sigma, mu, e):
    weights_record = []
    volatilities = []
    results = np.zeros((3,len(returns)))
    i = 0
    for portfolio_return in returns:
        A = np.block([[2*sigma, mu, e], [mu.T, 0, 0], [e.T, 0, 0]])
        b = np.zeros(n+2)
        b[n] = portfolio_return
        b[n+1] = 1
        w = np.linalg.solve(A, b)[:n]
        weights_record.append(w)
        portfolio_std_dev = np.sqrt( w.T @ sigma @ w )
        volatilities.append(portfolio_std_dev)
        results[0,i] = portfolio_std_dev
        results[1,i] = portfolio_return
        results[2,i] = (portfolio_return - risk_free_rate) / portfolio_std_dev
        i+=1
    return results, weights_record, volatilities

def display_simulated_ef_with_random(mean_returns, risk_free_rate, sigma, mu, e, df):
    results, weights, volatilities = efficient_portfolios(mean_returns,risk_free_rate, sigma, mu, e)
    max_sharpe_idx = np.argmax(results[2])
    sdp, rp = 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("-"*80)
    print("Maximum Sharpe Ratio Portfolio Allocation\n")
    print("Annualised Return:", round(rp,2))
    print("Annualised Volatility:", round(sdp,2))
    print("\n")
    print(max_sharpe_allocation)
    print("-"*80)
    print("Minimum Volatility Portfolio Allocation\n")
    print("Annualised Return:", round(rp_min,2))
    print("Annualised Volatility:", round(sdp_min,2))
    print("\n")
    print(min_vol_allocation)

    plt.figure(figsize=(10, 7))
    plt.scatter(results[0,:],results[1,:],c=results[2,:],cmap='YlGnBu', marker='o', s=10, alpha=0.3)
    plt.colorbar()
    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')
    plt.legend(labelspacing=0.8)

    return max_sharpe_allocation, min_vol_allocation

And it gives the following plot and portfolios:

--------------------------------------------------------------------------------
Maximum Sharpe Ratio Portfolio Allocation

Annualised Return: 0.63
Annualised Volatility: 0.23


               DD  ADBE  ATVI   APD    NVS      A   ADI    AVB    AYI   AAN  \
allocation -19.33  0.03 -0.32  29.3  12.65 -14.57  2.85 -25.28 -13.17  2.77   

            ...   SWKS    NOV  KMT   MDT   RIO   PSA   STE  POWI  VALE   TX  
allocation  ... -15.61 -10.08 -7.2 -3.16  7.57 -9.39  7.93  5.13  1.07  8.4  

[1 rows x 51 columns]
--------------------------------------------------------------------------------
Minimum Volatility Portfolio Allocation

Annualised Return: 0.03
Annualised Volatility: 0.13


             DD  ADBE  ATVI   APD   NVS     A   ADI   AVB   AYI   AAN  ...  \
allocation -0.6 -7.11  5.36  3.81  22.9 -3.69  7.37 -1.27 -1.13 -0.16  ...   

            SWKS   NOV   KMT   MDT   RIO    PSA   STE  POWI  VALE    TX  
allocation  -6.4 -0.25 -9.24  6.15  4.41  19.86 -1.31 -0.23 -2.99  6.05  

enter image description here

Let me know if I did anything wrong.

$\endgroup$
2
  • $\begingroup$ Are you trying to recover the means and returns from observed behavior or something? If the means and returns are given, you maximize the value of the portfolio and the weights are the solution: en.wikipedia.org/wiki/Modern_portfolio_theory $\endgroup$
    – user762914
    Jun 1, 2020 at 17:10
  • $\begingroup$ @Renard thanks for your comment. The returns are given indeed, I don't know about the means but I added a bit of context in the just created update section. I am trying to get the weights that would, for a given volatility $\sigma$ maximize the expected returns. I hope it makes sense? I am really not a mathematician $\endgroup$ Jun 1, 2020 at 17:27

1 Answer 1

1
+150
$\begingroup$

The following code generates points directly on the Pareto frontier. It uses the fact that the minimizer of a quadratic objective under linear constraints can be obtained by solving a linear system. The frontier changes shape if you uncomment the lines that add the risk-free asset.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# read data
stock_prices = pd.read_csv('Elegible stocks-Grid view.csv', index_col=0)
returns = stock_prices.pct_change()
mu = 252 * returns.mean().values
sigma = 252 * returns.cov().values
n = mu.shape[0]

# add risk free asset to mu/sigma
risk_free_rate = 0.0178
z = np.zeros((n,1))
#mu = np.block([mu, risk_free_rate])
#sigma = np.block([[sigma, z], [z.T, 0]])
#n = mu.shape[0]

# solve minimize w'∑w subject to μ'w = r, e'w=1 for varying r
mu = np.expand_dims(mu, axis=1)
e = np.ones((n,1))
returns = np.linspace(risk_free_rate, np.max(mu))
volatilities = []
for r in returns:
    A = np.block([[2*sigma, mu, e], [mu.T, 0, 0], [e.T, 0, 0]])
    b = np.zeros(n+2)
    b[n] = r
    b[n+1] = 1
    w = np.linalg.solve(A, b)[:n]

    volatility = np.sqrt( w.T @ sigma @ w )
    volatilities.append(volatility)

# plot results
plt.plot(volatilities, returns)
plt.xlabel('Volatility')
plt.ylabel('Return')
plt.show()
$\endgroup$
2
  • $\begingroup$ Thanks for your great answer. It is very clear and optimal, I will leave the question open for one more day if anybody has anything to add and will accept it. I had just two questions. 1. How would you if I needed w to be always be positive? 2. I didn't understood why the frontier changes shape when uncommenting the lines that add the risk-free asset. What do they stand for ? $\endgroup$ Jun 7, 2020 at 16:06
  • $\begingroup$ 1. then the solution cannot simply be obtained by solving a linear system. You would need to use a quadratic optimization solver; 2. I think the line is what is called best possible CAL on Wikipedia. $\endgroup$
    – LinAlg
    Jun 7, 2020 at 18:02

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .