import pandas as pd
import random
import math
from gurobipy import Model, quicksum, GRB, multidict
# from mypulp import Model, quicksum, GRB, multidict
from scipy.stats import norm
import numpy as np
import yfinance as yf
import warnings
import riskfolio as rp
import riskfolio.PlotFunctions as plf
Markowitzモデル
Markowitzモデルは，ポートフォリオ最適化問題の古典であり，以下のように定義される．
資産 $M$ 円を持つとき， $1$ 年後の期待資産価値が $\alpha M$ 円以上(ただし $\alpha >1$)という制約のもとで, 「リスク」を最小化することを考える． ただし，株を $i=1,2,\ldots,n $ で表し， 株 $i$ の $1$ 年後の価値は期待値 $r_i$, 分散 $\sigma^2_i$ の確率分布に したがうと仮定する．
株を $1, 2, \ldots, n$ とし， 各株 $i$ の $1$ 年後の価格を，確率変数 $R_i$ で表す． ただし，$R_i$ の期待値は $r_i$, 分散は $\sigma_i^2$ とする． 期待値をとる操作を $\mathbf{E}[\cdot]$ で書けば，次の関係が成り立っている． $$ \mathbf{E}[R_i] = r_i, \ \mathbf{E}[(R_i- r_i)^2] = \sigma_i^2 \ \ \forall i=1,2,\ldots,n $$
株 $i$ に投資する割合を $x_i$ とする．すると $x_1, x_2, \ldots, x_n$ は 非負で和が $1$ になっているはずである． $$ \sum_{i=1}^n x_i =1, \quad x_i\geq 0, \ \forall i=1,2,\ldots, n $$ この投資比率で投資したとき， $1$ 年後の財産価格は確率変数 $R_i$ を用いて $M \sum_{i=1}^n R_i x_i$ と書けるので，期待値が $ \alpha M$以上という制約は次のようになる． $$ M \sum_{i=1}^n r_i x_i \geq \alpha M $$
$M$ は両辺に現れるので消去でき，結局以下のようになる． $$ \sum_{i=1}^n r_i x_i \geq \alpha $$
最後に「リスク」とは何かを考えなければならない． Markowitz は「リスク」とは期待値からのずれと解釈し， 確率でいう「分散」をリスクと定義した．
$$ \begin{array}{ll} minimize &\displaystyle\sum_{i=1}^n \sigma_i^2 x_i^2 \\ s.t. & \displaystyle\sum_{i=1}^n r_i x_i \geq \alpha \\ & \displaystyle\sum_{i=1}^n x_i = 1 \\ & x_i \geq 0 \ \ \ \forall i=1, 2, \ldots, n \end{array} $$以下に定式化のコードと求解例を示す．
def markowitz(I, sigma, r, alpha):
"""markowitz -- simple markowitz model for portfolio optimization.
Parameters:
- I: set of items
- sigma[i]: standard deviation of item i
- r[i]: revenue of item i
- alpha: acceptance threshold
Returns a model, ready to be solved.
"""
model = Model("markowitz")
x = {}
for i in I:
x[i] = model.addVar(vtype="C", name="x(%s)" % i) # quantity of i to buy
model.update()
model.addConstr(quicksum(r[i] * x[i] for i in I) >= alpha)
model.addConstr(quicksum(x[i] for i in I) == 1)
model.setObjective(quicksum(sigma[i] ** 2 * x[i] * x[i] for i in I), GRB.MINIMIZE)
model.update()
model.__data = x
return model
I, sigma, r = multidict(
{1: [0.07, 1.01], 2: [0.09, 1.05], 3: [0.1, 1.08], 4: [0.2, 1.10], 5: [0.3, 1.20]}
)
alpha = 1.05
model = markowitz(I, sigma, r, alpha)
model.optimize()
x = model.__data
EPS = 1.0e-6
print("%5s\t%8s" % ("i", "x[i]"))
for i in I:
print("%5s\t%8g" % (i, x[i].X))
print("sum:", sum(x[i].X for i in I))
print("Obj:", model.ObjVal)
損をする確率を抑えるモデル
Markowitzのモデルは，リスクを分散で表現し， これを最小化する問題であった． しかし分散は，プラスの方向にもマイナスの方向にも等しく 作用するので，たくさん損することを避けるために， たくさん得することも避けるモデルとなっている． これは少し納得しづらい状況である．
そこで少し考え方を変える．もっと直接的に確率を扱い， 例えば $1$ 年後の資産価値が $0.95 M$ 円以下になる確率を $5\%$ 以下にしたいと考えてみよう． つまり，$95\%$ 以上の確率で 9割5分の資産が残る，という条件のもと， 期待資産価値を最大化する問題である．
上の問題を一般的に表すと次のようになる．
$1$ 年後の資産価値が $\alpha M$円以下になる「確率」を $\beta \%$ 以下に抑えながら， $1$ 年後の期待資産価値を最大化せよ.
先と同じように株式を $i=1,2,\ldots,n $ で表し， 株式 $i$ の一年後の価値が期待値 $r_i$, 分散 $\sigma_i$ の正規分布をすると仮定する． この確率変数を $R_i$ と書く． 各株式 $i$ に投資する割合を $x_i$ とすれば，やはり次を満たしている． $$ \sum_{i=1}^n x_i =1, \quad x_i\geq 0\ \ \ \forall i=1, 2, \ldots, n $$ 今回の場合，期待利益を最大化したいので，目的関数は $ \sum_{i=1}^n r_i x_i $ である．
難しいのは，「資産価値が $\alpha M$ 円以下になる確率を $\beta$ 以下に押さえる」という条件である．
確率変数 $R_i$ は平均 $r_i$, 分散 $\sigma_i^2$ の正規分布にしたがうと仮定したので， $R_i x_i$ は平均 $r_i x_i$, 分散 $\sigma_i^2 x_i^2$ の正規分布にしたがい， さらにその和 $ \sum_{i=1}^n R_i x_i $ は 平均 $\mu = \sum_{i=1}^n r_i x_i$, 分散 $\sigma^2 = \sum_{i=1}^n r_i^2 x_i^2$ の正規分布にしたがう． よって，新たな確率変数を $$ R = \frac{\sum_{i=1}^n R_i x_i -\mu}{\sigma} $$ で定義すれば $R$ は平均 $0$, 分散 $1$ の正規分布にしたがうことになる．
よって， $\Phi$ を正規分布の分布関数としたとき， $$ \begin{array}{l c l} Prob\left{ R \leq \frac{\alpha-\mu}{\sigma} \right} = \Phi\left(R \leq \frac{\alpha-\mu}{\sigma}\right) \leq \beta & \Leftrightarrow & \frac{\alpha - \mu}{\sigma} \leq \Phi^{-1}(\beta) \ & \Leftrightarrow &
- \Phi^{-1}(\beta)\sqrt{\sum_{i=1}^n\sigma_i^2 xi^2} \leq -\alpha + \sum{i=1}^n r_i x_i \end{array} $$ となる． $\beta<1/2$ のとき，$\Phi^{-1}(\beta)<0$ であり， 上の制約は2次錐制約となる．
結局，以下の最適化問題に定式化されることになる． $$ \begin{array}{ll} maximize & \rho \\ s.t. & \displaystyle\sum_{i=1}^n r_i x_i = \rho \\ & \displaystyle\sum_{i=1}^n x_i = 1 \\ & x_i \ge 0 \ \ \forall i=1,\ldots,n \\ & \sqrt{\displaystyle\sum_{i=1}^n \bar\sigma_i^2 x_i^2} \leq \displaystyle\frac{\alpha - \rho}{\Phi^{-1}(\beta)} \end{array} $$
これは2次錐最適化問題なので，Gurobiで解くことができる．
def p_portfolio(I, sigma, r, alpha, beta):
"""p_portfolio -- modified markowitz model for portfolio optimization.
Parameters:
- I: set of items
- sigma[i]: standard deviation of item i
- r[i]: revenue of item i
- alpha: acceptance threshold
- beta: desired confidence level
Returns a model, ready to be solved.
"""
model = Model("p_portfolio")
x = {}
for i in I:
x[i] = model.addVar(vtype="C", name="x(%s)" % i) # quantity of i to buy
rho = model.addVar(vtype="C", name="rho")
rhoaux = model.addVar(vtype="C", name="rhoaux")
model.update()
model.addConstr(rho == quicksum(r[i] * x[i] for i in I))
model.addConstr(quicksum(x[i] for i in I) == 1)
model.addConstr(rhoaux == (alpha - rho) / norm.ppf(beta))
model.addConstr(quicksum(sigma[i] ** 2 * x[i] * x[i] for i in I) <= rhoaux * rhoaux)
model.setObjective(rho, GRB.MAXIMIZE)
model.update()
model.__data = x
return model
I, sigma, r = multidict(
{1: [0.07, 1.01], 2: [0.09, 1.05], 3: [0.1, 1.08], 4: [0.2, 1.10], 5: [0.3, 1.20]}
)
alpha = 0.95
beta = 0.1
print("\n\n\nbeta:", beta, "phi inv:", phi_inv(beta))
model = p_portfolio(I, sigma, r, alpha, beta)
model.optimize()
x = model.__data
EPS = 1.0e-6
print("Investment:")
print("%5s\t%8s" % ("i", "x[i]"))
for i in I:
print("%5s\t%8g" % (i, x[i].X))
print("Obj:", model.ObjVal)
様々なモデルを解くためのパッケージ
以下のパッケージを使うと，様々なポートフォリオ最適化問題を解くことができる．
https://github.com/dcajasn/Riskfolio-Lib
データはYahoo Financeデータを yfinance パッケージを利用して入手している．ただし，現在はデータスクレイピングが禁止されたため，使えなくなっている． 使用する際には，自分でデータを準備する必要がある．
最適化はオープンソースの凸最適化パッケージ cvxopt を利用している．
最適化と主な引数
最適化は以下の形式で行う．
import riskfolio.Portfolio as pf
portofolio = pf.Portfolio()
portofolio.optimization()
主な引数は以下の通り．
model: 平均と共分散のためのモデル
- 'Classic': 過去の履歴（既定値）
- 'BL': Black Littermanモデル
- 'FM': リスクファクターモデル
rm: リスク尺度
'MV': 標準偏差 （既定値）
平均絶対偏差 (Mean Absolute Deviation: MAD)
- 'MSV': セミ標準偏差
- 'FLPM': First Lower Partial Moment (Omega Ratio)
- 'SLPM': Second Lower Partial Moment (Sortino Ratio)
- 'CVaR': 条件付きVaR
- 'EVaR': Entropic Value at Risk
- 'WR': Worst Realization (Minimax)
- 'MDD': Maximum Drawdown of uncompounded cumulative returns (Calmar Ratio)
- 'ADD': Average Drawdown of uncompounded cumulative returns
- 'CDaR': Conditional Drawdown at Risk of uncompounded cumulative returns
ただし， $$ \text{DaR}_{\alpha}(X) = \max_{j \in (0,T)} \left \{ \text{DD}(X,j) \in \mathbb{R}: F_{\text{DD}} \left ( \text{DD}(X,j) \right )< 1-\alpha \right \} $$
$$ \text{DD}(X,j) = \max_{t \in (0,j)} \left ( \sum_{i=0}^{t}X_{i} \right )- \sum_{i=0}^{j}X_{i} $$- 'EDaR': Entropic Drawdown at Risk of uncompounded cumulative returns
- 'UCI': Ulcer Index of uncompounded cumulative returns
obj: 目的関数
- 'MinRisk': リスク最小化
- 'Utility': 効用関数 $\mu w - l \phi_{i}(w)$ 最大化 （$\phi_{i}$ は上の13のリスク尺度から選択， $l$ はリスク回避率）
- 'Sharpe': リスク調整済みリターン（既定値）； (リターン - 無リスク金利)/$\phi_{i}$ で定義される．
- 'MaxRet': 期待リターン最大化
まず，データを読み込む．
warnings.filterwarnings("ignore")
pd.options.display.float_format = '{:.4%}'.format
# Date range
start = '2016-01-01'
end = '2019-12-30'
# Tickers of assets
assets = ['JCI', 'TGT', 'CMCSA', 'CPB', 'MO', 'APA', 'MMC', 'JPM',
'ZION', 'PSA', 'BAX', 'BMY', 'LUV', 'PCAR', 'TXT', 'TMO',
'DE', 'MSFT', 'HPQ', 'SEE', 'VZ', 'CNP', 'NI', 'T', 'BA']
assets.sort()
# Downloading data
data = yf.download(assets, start = start, end = end)
data = data.loc[:,('Adj Close', slice(None))]
data.columns = assets
Y = data[assets].pct_change().dropna()
Y.head()
ポートフォリオを求める．
port = rp.Portfolio(returns=Y)
method_mu='hist' # 過去の履歴から平均値を計算（ewma1 or ewma2にすると，指数平滑法）
method_cov='hist' # 過去の履歴から共分散を計算
port.assets_stats(method_mu=method_mu, method_cov=method_cov, d=0.94) #dは指数平滑法のパラメータ
# 最適ポートフォリオを求める
model='Classic' # Classic (historical), BL (Black Litterman) or FM (Factor Model)
rm = 'MV' # リスク尺度は標準偏差
obj = 'Sharpe' # 目的関数 MinRisk, MaxRet, Utility or Sharpe
hist = True # Use historical scenarios for risk measures that depend on scenarios
rf = 0 # Risk free rate
l = 0 # Risk aversion factor, only useful when obj is 'Utility'
w = port.optimization(model=model, rm=rm, obj=obj, rf=rf, l=l, hist=hist)
display(w.T)
結果を描画する．
ax = plf.plot_pie(w=w, title='Sharpe Mean Variance', others=0.05, nrow=25, cmap = "tab20",
height=6, width=10, ax=None)
points = 50 # Number of points of the frontier
frontier = port.efficient_frontier(model=model, rm=rm, points=points, rf=rf, hist=hist)
display(frontier.T.head())
label = 'Max Risk Adjusted Return Portfolio' # Title of point
mu = port.mu # Expected returns
cov = port.cov # Covariance matrix
returns = port.returns # Returns of the assets
ax = plf.plot_frontier(w_frontier=frontier, mu=mu, cov=cov, returns=returns, rm=rm,
rf=rf, alpha=0.05, cmap='viridis', w=w, label=label,
marker='*', s=16, c='r', height=6, width=10, ax=None)
ax = plf.plot_frontier_area(w_frontier=frontier, cmap="tab20", height=6, width=10, ax=None)
