現代ポートフォリオ理論と最大SRポートフォリオの分析

はじめに

本日は、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)@[e-mail](wt)=variance of the PF
def portfolio_std(wt, rt=None, cov =None):
    """Annualized PF standard deviation"""
    return np.sqrt((np.array(wt)@[e-mail](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