Portfolio optimization & backtesting
We evaluate, compare, and demonstrate different packages for performing portfolio optimization. There are several options available
- Optimization using scipy.optimize
- Optimization with cvxopt
- Optimiation with cvxpy
To compare the validity of our results, we will replicate the dataset and time window applied by DeMiguel et al. (2009) and its accompanying appendix.
We will replicate the Sharpe Ratios in in (Table 3) of DeMiguel et al. (2009). As long as the results are reasonably close, they will be deemed acceptable. This is because minor changes are made to the Ken French's dataset, and optimization is sensitive to changes in the input data.
Import modules¶
In [1]:
import os
import pandas as pd
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from pprint import pprint
from scipy.optimize import minimize
import pandas_datareader.data as web # module for reading datasets directly from the web
from pandas_datareader.famafrench import get_available_datasets
Reading data from Ken French's website¶
I print out all of the available datasets and find that there are 286 of them
In [2]:
datasets = get_available_datasets()
print('No. of datasets:{0}'.format(len(datasets)))
#datasets #comment out if you want to see all the datasets
I am looking for the 10 industry portfolio, and want to find which are the right keywords
In [3]:
df_10_industry = [dataset for dataset in datasets if '10' in dataset and 'Industry' in dataset]
print(df_10_industry)
Note that if you do not insert start or end dates, the default will obtain portfolios from approximately 2010 to the latest datapoints available.
In [4]:
ds_industry = web.DataReader(df_10_industry[0],'famafrench',start='1963-07-01',end='2004-11-01') # Taking [0] as extracting '10_Industry_Portfolios'
Obtaining data from the datareader returns a dict. Thus we want to see what is inside the dict.
In [5]:
print(type(ds_industry))
ds_industry.keys()
Out[5]:
We see that there is a DESCR. This tells us what the different datasets inside the dictionary refer to.
In [6]:
print(ds_industry['DESCR'])
In order to replicate DeMiguel et al. (2009), we will use dataset 0.
In [7]:
ds_industry[0].head()
Out[7]:
We perform the same process for the industry portfolios as for the FF risk factor portfolio dataset
In [8]:
df_5_factor = [dataset for dataset in datasets if '5' in dataset and 'Factor' in dataset]
print(df_5_factor)
ds_factors = web.DataReader(df_5_factor[0],'famafrench',start='1963-07-01',end='2004-11-01') # Taking [0] as extracting 1F-F-Research_Data_Factors_2x3')
print('\nKEYS\n{0}'.format(ds_factors.keys()))
print('DATASET DESCRIPTION \n {0}'.format(ds_factors['DESCR']))
ds_factors[0].head()
Out[8]:
In [9]:
dfAsset = ds_industry[0]/100
dfFactor = ds_factors[0]/100
Creating excess returns by subtracting the risk-free rate from the industry porfolio returns.
In [10]:
dfXsAsset = dfAsset.sub(dfFactor['RF'],axis=0)
dfXsAsset.head()
Out[10]:
Portfolio performance metrics¶
In [11]:
def sharpe(returns, rf=0, days=1):
volatility = returns.std() * np.sqrt(days)
sharpe_ratio = (returns.mean() - rf) / volatility
return sharpe_ratio
def information_ratio(returns, benchmark_returns, days=252):
return_difference = returns - benchmark_returns
volatility = return_difference.std() * np.sqrt(days)
information_ratio = return_difference.mean() / volatility
return information_ratio
Equal-weighted portfolio (1/N) (re-balancing)¶
In [12]:
M = 120
pfRet=[]
pfStratRet = dict()
pfSharpeR = dict()
for date,assRet in dfXsAsset[M:].iterrows():
N = len(assRet)
tgtWtg = 1/N*np.ones(N)
assRet = np.array(assRet.values)
pfRet.append(float(tgtWtg@assRet))
dfStratRet = pd.DataFrame(pfRet,index=dfXsAsset[M:].index,columns=['ew-rebal'])
In [13]:
dfStratRet.head()
Out[13]:
Equal-weighted portfolio (1/N) (buy & hold)
In [ ]:
for date,assRet in dfXsAsset[M:].iterrows():
N = len(assRet)
tgtWtg = 1/N*np.ones(N)
assRet = np.array(assRet.values)
pfRet.append(float(tgtWtg@assRet))
dfStratRet = pd.DataFrame(pfRet,index=dfXsAsset[M:].index,columns=['ew-rebal'])
In [38]:
dfStratRet.head()
Out[38]:
Mean-variance portfolio, Out-of-sample, Matrix Calculus¶
In [14]:
mean_OOS = dfXsAsset.rolling(M,min_periods=M,win_type=None).mean().dropna()
cov_OOS = dfXsAsset.rolling(M,min_periods=M,win_type=None).cov().dropna()
Target weights¶
In [15]:
pfTgtWts = []
for date, cov_df in cov_OOS.groupby(level=0):
N = cov_df.shape[0]
cov = cov_df.values
mean = mean_OOS.loc[date,:].values
inv_cov = np.linalg.pinv(cov)
tgtWts_num = inv_cov@mean
tgtWts_denom = np.ones((1,N))@inv_cov@mean
pfTgtWts.append(tgtWts_num/tgtWts_denom)
pfTgtWts = np.array(pfTgtWts)
pfTgtWts = pfTgtWts[:-1,:]
dfTgtWts = pd.DataFrame(pfTgtWts,index=dfXsAsset[M:].index,columns=list(cov_df))
dfTgtWts['Strategy'] = 'mv_oos'
dfTgtWts.set_index('Strategy',append=True,inplace=True)
Portfolio returns¶
In [16]:
pfRet = []
for date,assRet in dfXsAsset[M:].iterrows():
N = len(assRet)
tgtWtg = dfTgtWts.loc[date,:].values
pfRet.append(float(tgtWtg@assRet.values))
dfStratRet['mv_oos'] = pfRet
Mean-variance portfolio, Out-of-sample, Optimization (scipy package)¶
In [17]:
mean_OOS = dfXsAsset.rolling(M,min_periods=M,win_type=None).mean().dropna()
cov_OOS = dfXsAsset.rolling(M,min_periods=M,win_type=None).cov().dropna()
Target weights¶
In [18]:
# Calculates portfolio mean return
def port_mean(W, R):
return sum(R * W)
# Calculates portfolio variance of returns
def port_var(W, C):
return np.dot(np.dot(W, C), W)
def port_mean_var(W, R, C):
return port_mean(W, R), port_var(W, C)
def mvo_objFunc(W, R, C, rf):
mean, var = port_mean_var(W, R, C)
util = (mean - rf) / np.sqrt(var)
return 1 / util
def opt_weights(R, C, rf):
n = len(R)
W = np.ones([n]) / n
b_ = [(0, 1.0) for i in range(n)]
c_ = ({'type': 'eq', 'fun': lambda W: np.sum(W)-1.0})
optimized = minimize(mvo_objFunc, W, (R, C, rf), method='SLSQP', constraints=c_, bounds=b_)
return optimized.x
In [19]:
pfTgtWts = []
for date, cov_df in cov_OOS.groupby(level=0):
N = cov_df.shape[0]
cov = cov_df.values
mean = mean_OOS.loc[date,:].values
pfTgtWts.append(opt_weights(mean,cov,0))
In [20]:
pfTgtWts = np.array(pfTgtWts)
pfTgtWts[pfTgtWts<10**-6] = 0
pfTgtWts = pfTgtWts[:-1,:]
In [21]:
dfTgtWts = pd.DataFrame(pfTgtWts,index=dfXsAsset[M:].index,columns=list(cov_df))
dfTgtWts['Strategy'] = 'mv_oos_opt'
dfTgtWts.set_index('Strategy',append=True,inplace=True)
Portfolio returns¶
In [22]:
pfRet = []
for date,assRet in dfXsAsset[M:].iterrows():
N = len(assRet)
tgtWtg = dfTgtWts.loc[date,:].values
pfRet.append(float(tgtWtg@assRet.values))
dfStratRet['mv_oos_opt'] = pfRet
Mean-variance portfolio, Out-of-sample, Optimization (cvopt package)¶
In [23]:
mean_OOS = dfXsAsset.rolling(M,min_periods=M,win_type=None).mean().dropna()
cov_OOS = dfXsAsset.rolling(M,min_periods=M,win_type=None).cov().dropna()
Target weights¶
In [24]:
import cvxopt as opt
from cvxopt import blas, solvers
def optimal_portfolio(meanRet, cov):
n = len(meanRet)
N = 100 #number of points on efficient frontier
mus = [10**(5.0 * t/N - 1.0) for t in range(N)] # mean returns for efficient frontier
# Convert to cvxopt matrices
S = opt.matrix(cov)
pbar = opt.matrix(meanRet)
# Create constraint matrices
G = -opt.matrix(np.eye(n)) # negative n x n identity matrix
h = opt.matrix(0.0, (n ,1))
A = opt.matrix(1.0, (1, n))
b = opt.matrix(1.0)
# Calculate efficient frontier weights using quadratic programming
portfolios = [solvers.qp(mu*S, -pbar, G, h, A, b)['x'] for mu in mus]
## CALCULATE RISKS AND RETURNS FOR FRONTIER
returns = [blas.dot(pbar, x) for x in portfolios]
risks = [np.sqrt(blas.dot(x, S*x)) for x in portfolios]
## CALCULATE THE 2ND DEGREE POLYNOMIAL OF THE FRONTIER CURVE
m1 = np.polyfit(returns, risks, 2)
x1 = np.sqrt(m1[2] / m1[0])
# CALCULATE THE OPTIMAL PORTFOLIO
solvers.options['show_progress']= False
wt = solvers.qp(opt.matrix(x1 * S), -pbar, G, h, A, b)['x']
return wt
In [25]:
pfTgtWts = []
for date, cov_df in cov_OOS.groupby(level=0):
N = cov_df.shape[0]
cov = cov_df.values
mean = mean_OOS.loc[date,:].values
wtg = optimal_portfolio(mean,cov)
wtg = list(wtg)
wtg = np.asarray(wtg)
pfTgtWts.append(wtg)
In [26]:
pfTgtWts = np.array(pfTgtWts)
pfTgtWts[pfTgtWts<10**-6] = 0
pfTgtWts = pfTgtWts[:-1,:]
In [27]:
dfTgtWts = pd.DataFrame(pfTgtWts,index=dfXsAsset[M:].index,columns=list(cov_df))
dfTgtWts['Strategy'] = 'mv_oos_cvxopt'
dfTgtWts.set_index('Strategy',append=True,inplace=True)
Portfolio returns¶
In [28]:
pfRet = []
for date,assRet in dfXsAsset[M:].iterrows():
N = len(assRet)
tgtWtg = dfTgtWts.loc[date,:].values
pfRet.append(float(tgtWtg@assRet.values))
dfStratRet['mv_oos_cvxopt'] = pfRet
Portfolio results¶
The results below show that our ew (1/N) portfolio is very similar to Table 3 results, and so is the mv_oos_cvxopt. Only the implementation in scipy.optimize.minimize is quite different and we will need to investigate further
In [29]:
dfStratRet.apply(sharpe,axis=0)
Out[29]:
In [30]:
dfStratRet.describe()
Out[30]:
In [31]:
dfStratRet.head()
Out[31]:
In [37]:
dfStratRet['mv_oos'].mean()
Out[37]:
In [32]:
dfStratRet.plot(title='Return series')
Out[32]:
In [33]:
((1*dfStratRet+1).cumprod()-1).plot(title='Cumulative Wealth')
Out[33]:
In [ ]:
Comments
Comments powered by Disqus