Markovitz Portfolio Optimization in a simple example¶
The example follows closely http://www.pythonforfinance.net/2017/01/21/investment-portfolio-optimisation-with-python/ .
In [1]:
import matplotlib.pyplot as plt
import pandas_datareader.data as web
import numpy as np
import pandas as pd
start = pd.datetime(2013,1,1,0,0,0,0)
In [4]:
import numpy as np
import pandas as pd
import pandas_datareader.data as web
import matplotlib.pyplot as plt
#list of stocks in portfolio
stocks = ['AAPL','AMZN','MSFT','GOOG']
#download daily price data for each of the stocks in the portfolio
data = web.DataReader(stocks,data_source='yahoo',start=start)['Adj Close']
#data can be stored by pd.Dataframe.to_csv(data,'data.csv')
#data = pd.DataFrame.from_csv('data.csv')
#convert daily stock prices into daily returns
returns = -data.pct_change()
#calculate mean daily return and covariance of daily returns
mean_daily_returns = returns.mean()
cov_matrix = returns.cov()
#set number of runs of random portfolio weights
num_portfolios = 25000
#set up array to hold results
#We have increased the size of the array to hold the weight values for each stock
results = np.zeros((4+len(stocks)-1,num_portfolios))
for i in range(num_portfolios):
#select random weights for portfolio holdings
weights = np.array(np.random.random(4))
#rebalance weights to sum to 1
weights /= np.sum(weights)
#calculate portfolio return and volatility
portfolio_return = np.sum(mean_daily_returns * weights) * 252
portfolio_std_dev = np.sqrt(np.dot(weights.T,np.dot(cov_matrix, weights))) * np.sqrt(252)
#store results in results array
results[0,i] = portfolio_return
results[1,i] = portfolio_std_dev
#store Sharpe Ratio (return / volatility) - risk free rate element excluded for simplicity
results[2,i] = results[0,i] / results[1,i]
#iterate through the weight vector and add data to results array
for j in range(len(weights)):
results[j+3,i] = weights[j]
#convert results array to Pandas DataFrame
results_frame = pd.DataFrame(results.T,columns=['ret','stdev','sharpe',stocks[0],stocks[1],stocks[2],stocks[3]])
#locate position of portfolio with highest Sharpe Ratio
max_sharpe_port = results_frame.iloc[results_frame['sharpe'].idxmax()]
#locate positon of portfolio with minimum standard deviation
min_vol_port = results_frame.iloc[results_frame['stdev'].idxmin()]
%matplotlib inline
#create scatter plot coloured by Sharpe Ratio
plt.scatter(results_frame.stdev,results_frame.ret,c=results_frame.sharpe,cmap='RdYlBu',s=1)
plt.xlabel('Volatility')
plt.ylabel('Returns')
plt.colorbar()
#plot red star to highlight position of portfolio with highest Sharpe Ratio
plt.scatter(max_sharpe_port[1],max_sharpe_port[0],marker=(5,1,0),color='r',s=100)
#plot green star to highlight position of minimum variance portfolio
plt.scatter(min_vol_port[1],min_vol_port[0],marker=(5,1,0),color='g',s=100)
Out[4]:
In [5]:
data.plot()
Out[5]:
In [6]:
returns
Out[6]:
In [7]:
np.sqrt(cov_matrix*252)
Out[7]:
In [8]:
mean_daily_returns
Out[8]:
In [9]:
min_vol_port
Out[9]:
In [10]:
max_sharpe_port
Out[10]:
In [11]:
#we calculate the value of a portfolio with given weights
weights_max_sharpe = [max_sharpe_port[3],max_sharpe_port[4],max_sharpe_port[5],max_sharpe_port[6]]
weights_min_vol = [min_vol_port[3], min_vol_port[4], min_vol_port[5], min_vol_port[6]]
def portfolio_value_plot(weights,data):
portfolio_value = [1]
returns = -data.pct_change()
for i in range(1,len(data.AMZN)):
portfolio_value = portfolio_value + [portfolio_value[-1]+portfolio_value[-1]*(returns.AAPL[i]*weights[0]+returns.AMZN[i]*weights[1]+
returns.GOOG[i]*weights[2]+returns.MSFT[i]*weights[3])]
plt.plot(portfolio_value,linewidth=0.2)
def portfolio_value_plot_max(weights,data):
portfolio_value = [1]
returns = -data.pct_change()
for i in range(1,len(data.AMZN)):
portfolio_value = portfolio_value + [portfolio_value[-1]+portfolio_value[-1]*(returns.AAPL[i]*weights[0]+returns.AMZN[i]*weights[1]+
returns.GOOG[i]*weights[2]+returns.MSFT[i]*weights[3])]
plt.plot(portfolio_value,'r')
portfolio_value_plot_max(weights_max_sharpe,data)
portfolio_value_plot_max(weights_min_vol,data)
for i in range(10):
randomweights = np.array(np.random.random(4))
#rebalance weights to sum to 1
randomweights /= np.sum(randomweights)
portfolio_value_plot(randomweights,data)
plt.show()
