はじめに
本日は、Stefan Jasen氏によって書かれた"Machine Learning for Algorithmic trading, 2nd edition"のChapter 5, Mean-variance optimizationについて勉強した内容をまとめます。Pythonを用いたコーディングに際して、カナヲ定量分析さんの記事を参考にさせて頂きました。Stefan氏の解説だけでは理解できない部分を補うことができました。ありがとうございます。
現代ポートフォリオ理論の考え方
- 同じ期待リターンならvolが小さいPFを選択することが、もっとも好ましい。
- 各期待リターンにおいてvolが最も小さいPFの集合を効率的フロンティアと呼ぶ。
コードの概要
- yahoofinanceからSP500にリストされている30の株（クリーニングの過程で27に絞られる）を、Fredからrisk-free rateを取得。
- ディリクレ分布より重みを1万通りサンプリングし、それらを用いて運用した結果をプロット。また、最大SRPF・最小分散PFをプロット。
- scipy.optimizeからminimizeを用いてSRが最大化する重みを推定。また、その重みを用いた運用結果を計算。
#Get data from yfinance using pandas_DataReader
import numpy as np
import pandas as pd
import pandas_datareader as pdr
import datetime
import matplotlib.pyplot as plt
import seaborn as sns
from numpy.linalg import inv
from numpy.random import random, uniform, dirichlet, choice
from scipy.optimize import minimize
tickers = ['MMM','ABT','ACN','AFL','ARE','GOOGL','AMZN','AME','APA','AAPL'
,'BKR','CPB','KMX','CCL','DHR','CVS','ECL','EA','EIX','FOXA','BEN'
,'FCX','HCA','IP','KEY','LNC','MTCH','MS','NLSN','PVH']
df = pdr.DataReader(tickers, data_source ='yahoo', start ='2010-01-01', end = '2021-01-01');
#delete any rows with NaN
df = df.dropna(how = 'any',axis = 1);
#clean the df
df = df['Adj Close'];
#compute weekly returns
weekly_returns = df.resample('W').last().pct_change().dropna(axis=0, how ='all');
#interpret stocks as the indices(['MMM',...,'PVH'])
stocks = weekly_returns.columns
#the number of data points and stocks
n_obs, n_assets = weekly_returns.shape
#the number of data points and stocks
n_obs, n_assets = weekly_returns.shape
#the number of simulation
NUM_PF = 10000
#set the range of portfolio weights
from numpy.random import random, uniform, dirichlet, choice
x0 = uniform(0,1,n_assets)
#x0 = xo / np.sum(np.abs(x0))
x0 /=np.sum(np.abs(x0))
#一年に何回リターンが計測されているか平均を出す
periods_per_year = round(weekly_returns.resample('A').size().mean())
periods_per_year
#compute average return, covariance, inverse covariance matrix
mean_returns = weekly_returns.mean()
cov_matrix = weekly_returns.cov()
inverse_cov_matrix = pd.DataFrame(inv(cov_matrix), index = stocks, columns = stocks)
#compute risk-free rate using Fred
df2 = pdr.DataReader('^TNX', data_source ='yahoo', start ='2010-01-01', end = '2021-01-01');
df2 = df2['Adj Close']
#リターン一つ分に対応するリスクフリーリターン
treasury_10yr_monthly = (df2.resample('M').last().div(periods_per_year).div(100).squeeze())
rf_rate = treasury_10yr_monthly.mean()
#simulate random portfolio
def simulate_portfolios(mean_ret, cov, rf_rate = rf_rate, short=True):
alpha = np.full(shape=n_assets, fill_value=.05) #0.05がn-assets分含まれたndarrayを作成
weights = dirichlet(alpha=alpha, size=NUM_PF) #各アセットの係数の組み合わせをNUM_PF分作成
if short:
weights *= choice([-1,1],size=weights.shape)#-1
returns = weights @ mean_ret.values +1 #重みと平均リターン行列の積
returns = returns ** periods_per_year -1#年率リターン
std = (weights @ weekly_returns.T).std(1)#週ごとのvolを計算
std *=np.sqrt(periods_per_year)#週率リスクから年率リスクへ変換
sharpe = (returns - rf_rate)/std
return pd.DataFrame({'Annualized Standard Deviation':std,
'Annualized Returns':returns,
'Sharpe Ratio': sharpe}), weights;
#simul_perfにはdfを、simul_wtにはalpha=0.05のディリクレ分布から取得した係数を格納
simul_perf, simul_wt = simulate_portfolios(mean_returns, cov_matrix, short=False);
df = pd.DataFrame(simul_wt);
#axesおよびx-axis,y-axisを設定する
axes = simul_perf.plot(kind='scatter', x="Annualized Standard Deviation",
y="Annualized Returns",c = 2, cmap ='Blues',alpha=0.5, figsize=(14,9), colorbar=True,
title='{} Simulated Portfolios'.format(NUM_PF))
#simul_perfから最大sharpe ratioをマークしたインデックスをilocおよびidxmax()を用いて取得
max_sharpe_idx = simul_perf.iloc[:,2].idxmax()
#iloc[max_sharpe_idx,:2]は二列目を含まないことに注意
sd, r = simul_perf.iloc[max_sharpe_idx, :2].values
print('Max Sharpe:{:.2%},{:.2%}'.format(sd,r))
⇒Max Sharpe:22.77%,34.30%
axes.scatter(x=sd,y=r, marker ='*', color = 'darkblue', s = 500, label = 'Max.Sharpe Ratio')
#simul_perfから最小リスクをマークしたインデックスをilocおよびidxmax()を用いて取得
min_vol_idx = simul_perf.iloc[:,0].idxmin()
sd,r = simul_perf.iloc[min_vol_idx,:2].values
axes.scatter(x=sd, y=r, marker ='*', color ='green', s = 500, label='Min Volatility')
plt.legend()
#最大シャープPFのパフォーマンスを評価するためにメソッドを定義する
#np.array(wt)@cov_values@np.array(wt)=variance of the PF
def portfolio_std(wt, rt=None, cov =None):
"""Annualized PF standard deviation"""
return np.sqrt((np.array(wt)@cov.values@np.array(wt))*periods_per_year)
#(np.array(wt)@np.array(rt))+1=PFのリターン（これを投資額と掛ければ元金を含めたリターンが出せる）
def portfolio_returns(wt, rt=None):
"""Annualized PF returns"""
return ((np.array(wt)@np.array(rt))+1)
def portfolio_performance(wt,rt,cov):
"""Annualized PF returns & standard deviation"""
r = portfolio_returns(wt, rt=rt)
sd = portfolio_std(wt,cov=cov)
return r,sd
#最大シャープPFを計算する関数を定義
def neg_sharpe_ratio(weights, mean_ret,cov):
r,sd = portfolio_performance(weights, mean_ret, cov)
return -(r-rf_rate)/sd
weights_constraint = {'type':'eq', 'fun':lambda x : np.sum(np.abs(x))-1}
#fun=最小化する関数、x0 = 最適化単探索の際に動かすパラメータ、args = パラメータ以外の入力値
#constraints bounds = パラメタへの制約
#options = 許容誤差と最大イテレーション回数
def max_sharpe_ratio(mean_ret, cov, short=False):
return minimize(fun=neg_sharpe_ratio,
x0=x0,
args=(mean_ret,cov),
method='SLSQP',
bounds=((-1 if short else 0,1),)*n_assets,
constraints=weights_constraint,
options={'ftol':1e-10,'maxiter':1e4})
max_sharpe_pf = max_sharpe_ratio(mean_returns, cov_matrix, short=False)
#max_sharpe_pf.x=最適化された重みの配列
max_sharpe_perf = portfolio_performance(max_sharpe_pf.x, mean_returns, cov_matrix)
r, sd = max_sharpe_perf
#indexとdataのみ
pd.Series({'ret': r, 'sd': sd, 'sr': (r-rf_rate)/sd})
ret 1.002794
sd 0.138476
sr 7.238293
dtype: float64
