Portfolio Optimization¶
“Modern Portfolio Theory (MPT), a hypothesis put forth by Harry Markowitz in his paper “Portfolio Selection,” (published in 1952 by the Journal of Finance) is an investment theory based on the idea that risk-averse investors can construct portfolios to optimize or maximize expected return based on a given level of market risk, emphasizing that risk is an inherent part of higher reward. It is one of the most important and influential economic theories dealing with finance and investment.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
Get Daily Returns from CSV Files¶
aapl = pd.read_csv('AAPL_CLOSE',index_col='Date', parse_dates=True)
cisco = pd.read_csv('CISCO_CLOSE',index_col='Date', parse_dates=True)
ibm = pd.read_csv('IBM_CLOSE',index_col='Date', parse_dates=True)
amzn = pd.read_csv('AMZN_CLOSE',index_col='Date', parse_dates=True)
Consolidate Columns into One DataFrame¶
Use concat() to pull columns axis=1 together in one DataFrame
stocks = pd.concat([aapl, cisco, ibm, amzn],axis=1)
Rename the columns by using member variable .columns
stocks.columns = ['aapl', 'cisco', 'ibm', 'amzn']
stocks.head()
Get Average Daily Return¶
Use .pct_change().mean to get the mean daily return
mean_daily_ret = stocks.pct_change(1).mean()
mean_daily_ret
Get Correlation Between the Returns¶
Use .pct_change().corr() to get the correlation
stocks.pct_change(1).corr()
Random Optimization¶
stocks.head()
stock_normed = stocks/stocks.iloc[0]
stock_normed.plot()
stock_daily_ret = stocks.pct_change(1)
stock_daily_ret.head()
Logarithmic Returns vs Arithmetic Returns¶
We will now switch over to using log returns instead of arithmetic returns, for many of our use cases they are almost the same, but most technical analyses require detrending/normalizing the time series and using log returns is a nice way to do that. Log returns are convenient to work with in many of the algorithms we will encounter.
For a full analysis of why we use log returns, check this great article.
Arithmetic Returns
stocks.pct_change(1).head()
Switch to Log Returns¶
log_ret = np.log(stocks/stocks.shift(1))
log_ret.head()
Plot Histograms¶
log_ret.hist(bins=100,figsize=(12,6));
plt.tight_layout()
Get the Mean Return¶
log_ret.mean()
log_ret.describe().transpose()
log_ret.mean() * 252
Get Pairwise Covariance of Columns¶
log_ret.cov()
log_ret.cov()*252 # multiply by days
Single Run for Some Random Allocation¶
# Set seed (optional)
# Get the same random num every time
np.random.seed(101)
# Stock Columns
print('Stocks')
print(stocks.columns)
print('\n')
# Create Random Weights
print('Creating Random Weights')
weights = np.array(np.random.random(4))
print(weights)
print('\n')
# Rebalance Weights
print('Rebalance to sum to 1.0')
weights = weights / np.sum(weights)
print(weights)
print('\n')
# Expected Return
print('Expected Portfolio Return')
exp_ret = np.sum(log_ret.mean() * weights) *252
print(exp_ret)
print('\n')
# Risk Free Return
# zero
# Expected Variance
# Linear Algebra is used here to make the code run faster
print('Expected Volatility')
exp_vol = np.sqrt(np.dot(weights.T, np.dot(log_ret.cov() * 252, weights)))
print(exp_vol)
print('\n')
# Sharpe Ratio
SR = exp_ret/exp_vol
print('Sharpe Ratio')
print(SR)
Simulating Thousands of Possible Allocations¶
Great! Now we can just run this many times over! First, we will remove the print() statements
# np.random.seed(101)
num_ports = 15000
all_weights = np.zeros((num_ports, len(stocks.columns)))
ret_arr = np.zeros(num_ports)
vol_arr = np.zeros(num_ports)
sharpe_arr = np.zeros(num_ports)
for ind in range(num_ports):
# Create Random Weights
weights = np.array(np.random.random(4))
# Rebalance Weights
weights = weights / np.sum(weights)
# Save Weights
all_weights[ind,:] = weights
# Expected Return
ret_arr[ind] = np.sum((log_ret.mean() * weights) *252)
# Expected Variance
vol_arr[ind] = np.sqrt(np.dot(weights.T, np.dot(log_ret.cov() * 252, weights)))
# Sharpe Ratio
sharpe_arr[ind] = ret_arr[ind]/vol_arr[ind]
Get Max. Sharpe Ratio
sharpe_arr.max()
Get Location of Max. Sharpe Ratio
sharpe_arr.argmax()
Get the Optimal Allocations
all_weights[1419,:]
Get the Return & Volatiolity at Max: used in the plot below
max_sr_ret = ret_arr[1419]
max_sr_vol = vol_arr[1419]
Plotting the data¶
plt.figure(figsize=(12,8))
plt.scatter(vol_arr, ret_arr, c=sharpe_arr, cmap='plasma')
plt.colorbar(label='Sharpe Ratio')
plt.xlabel('Volatility')
plt.ylabel('Return')
# Add red dot for max SR
plt.scatter(max_sr_vol, max_sr_ret, c='red', s=50, edgecolors='black')
Mathematical Optimization¶
There are much better ways to find good allocation weights than just guess and check! We can use optimization functions to find the ideal weights mathematically!
Functionalize Return and SR operations¶
def get_ret_vol_sr(weights):
"""
Takes in weights, returns array or return,volatility, sharpe ratio
"""
weights = np.array(weights)
ret = np.sum(log_ret.mean() * weights) * 252
vol = np.sqrt(np.dot(weights.T, np.dot(log_ret.cov() * 252, weights)))
sr = ret/vol
return np.array([ret,vol,sr])
Import scipy¶
from scipy.optimize import minimize
To fully understand all the parameters, check out: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html
Using help() to Get Function Help¶
help(minimize)
Optimization works as a minimization function, since we actually want to maximize the Sharpe Ratio, we will need to turn it negative so we can minimize the negative sharpe (same as maximizing the postive sharpe)
Define Help Functions¶
def neg_sharpe(weights):
return get_ret_vol_sr(weights)[2] * -1
# Contraints
def check_sum(weights):
'''
Returns 0 if sum of weights is 1.0
'''
return np.sum(weights) - 1
# By convention of minimize function it should be a function that returns zero for conditions
cons = ({'type':'eq','fun': check_sum})
# 0-1 bounds for each weight
bounds = ((0, 1), (0, 1), (0, 1), (0, 1))
# Initial Guess (equal distribution)
init_guess = [0.25,0.25,0.25,0.25]
Run Minimize() Function¶
# Sequential Least SQuares Programming (SLSQP).
opt_results = minimize(neg_sharpe, init_guess, method='SLSQP', bounds=bounds, constraints=cons)
opt_results
Get the Optimized Allocations¶
opt_results.x
Get the Optimal Results for Return, Volatility & Sharpe Ratio¶
get_ret_vol_sr(opt_results.x)
All Optimal Portfolios (Efficient Frontier)¶
The efficient frontier is the set of optimal portfolios that offers the highest expected return for a defined level of risk or the lowest risk for a given level of expected return. Portfolios that lie below the efficient frontier are sub-optimal, because they do not provide enough return for the level of risk. Portfolios that cluster to the right of the efficient frontier are also sub-optimal, because they have a higher level of risk for the defined rate of return.
Efficient Frontier http://www.investopedia.com/terms/e/efficientfrontier
From the plot above, we notice our returns go from 0 to somewhere along 0.3. Use linspace() to create a linspace number of points to calculate x on
# Change 100 to a lower number for slower computers!
frontier_y = np.linspace(0, 0.3, 100)
def minimize_volatility(weights):
return get_ret_vol_sr(weights)[1]
frontier_volatility = []
for possible_return in frontier_y:
# function for return
cons = ({'type':'eq','fun': check_sum},
{'type':'eq','fun': lambda w: get_ret_vol_sr(w)[0] - possible_return})
result = minimize(minimize_volatility, init_guess, method='SLSQP', bounds=bounds, constraints=cons)
frontier_volatility.append(result['fun'])
Plot Efficient Frontier Line¶
plt.figure(figsize=(12,8))
plt.scatter(vol_arr, ret_arr, c=sharpe_arr, cmap='plasma')
plt.colorbar(label='Sharpe Ratio')
plt.xlabel('Volatility')
plt.ylabel('Return')
# Add frontier line
plt.plot(frontier_volatility, frontier_y, 'g--', linewidth=3)
