はじめに

本記事では，リクルート社で開発された最適化フレームワークである「Codable Model Optimizer」を使ってポートフォリオ最適化問題を解いてみました．

最適化問題とは，膨大な選択肢からより良い選択を見つけ出す問題を指し，様々な解法があります．解法としては，厳密解法（主に厳密な最適解を見つけるのに使われる）やメタヒューリスティクス（厳密な最適解ではないが，大規模な問題において良い解を見つけるのに使われる）などがありますが，今回紹介する「Codable Model Optimizer」は数理最適化のプロではなくてもITエンジニアが気軽に利用できる最適化フレームワークを目指し開発されたようです．

現代には様々な社会課題に最適化問題が内在しており，最適化問題を解くことはSociety5.0やDX化が進む現代においてさらに重要になってくることが予想されます．今回はそんな最適化問題の中からポートフォリオ最適化問題に着目し「Codable Model Optimizer」を使って最適化を実行してみたので紹介いたします．

「Codable Model Optimizer」の特徴や使い方の詳細はこちらの記事で紹介されています．
参考記事：https://blog.recruit.co.jp/data/articles/codable_optimizer/
また，ポートフォリオ最適化問題を取り扱うにあたり，以下の記事を参考にさせていただきました．
参考記事
-https://qiita.com/yumaloop/items/d709cc9b43f18df70382
-https://qiita.com/well_living/items/b323a363e5442a15db6d

ポートフォリオ最適化問題とは

ポートフォリオ最適化とは複数の金融資産に対して適当な投資比率をそれぞれ決定することを指し，ポートフォリオ理論に基づき最適化が行われます．ポートフォリオ理論とは、金融資産の収益率に着眼点を置き，リターンとリスクを基にポートフォリオを決定する理論になりますが，多くの投資家はリスクを最小化しながらリターンを得るという考えを共通としてポートフォリオを決定します．今回は，このような投資家の合理性を反映した最も基本的なポートフォリオ最適化として知られている，Markowitzの平均分散モデル(Mean-Variance Model)を採用し最適化を行っていきます．

H. M. Markowitzの平均・分散モデル

Markowitzの平均分散モデルでは，「ポートフォリオの期待収益率が一定値（目標収益率）以上となる」という制約条件の下で，「ポートフォリオの分散（リスク）を最小化する」最適化問題を考えます．

一般に，$n$コの金融資産で構成されるポートフォリオの場合，ポートフォリオの分散は$n$コの資産間の共分散行列の二次形式となるので，この最適化問題は二次計画問題(Quadratic Programming, QP)のクラスとなり，次のように定式化されます．


\begin{align}
\underset{\bf x}{\rm minimize} ~~~ &{\bf x}^T \Sigma {\bf x} \\\
{\rm subject~to} ~~~ &{\bf r}^T {\bf x} = \sum_{i=1}^{n} r_i x_i \geq r_e \\\
&{\|{\bf x}\|}_{1} = \sum_{i=1}^{n} x_i = 1 \\\
&x_i \geq 0 ~~ (i = 1, \cdots, n)
\end{align}

$\Sigma \in \mathbb{R}^{n \times n}$ ー $n$コの金融資産の共分散行列
${\bf x} \in \mathbb{R}^{n}$ ー $n$コの金融資産の投資比率ベクトル
${\bf r} \in \mathbb{R}^{n}$ ー $n$コの金融資産の期待収益率ベクトル
$x_i \in \mathbb{R}$ ー金融資産$i$の投資比率
${r}_i \in \mathbb{R}$ ー金融資産$i$の期待収益率
$r_e \in \mathbb{R}$ ー投資家の要求期待収益率（目標収益率）

ここで，ポートフォリオの期待収益率は様々な方法が存在しますが，今回は時系列データの平均値を用いています．
その上で，各制約式についての補足を以下に示します．
制約式1：ポートフォリオの期待収益率が一定値($=r_e$)以上となることを要請．
制約式2：ポートフォリオの定義から自明．
制約式3：ポートフォリオの定義から自明．資産の空売りを許す場合，制約式3は除くこともある．

それでは，Pythonを用いて，実際に実装していきます．

Pythonで実装

準備

まずは最適化フレームワークである「Codable Model Optimizer」や株価データの取得に必要な「pandas_datareader」を環境にインストールします．

# Install codableopt
pip install codableopt
# Install pandas-datareader
pip install pandas-datareader

続いて，今回用いるライブラリをインポートしていきます．

import datetime
import numpy as np
import pandas as pd
import pandas_datareader.stooq as stooq
import matplotlib.pyplot as plt
from codableopt import *

株価データの準備

ポートフォリオ最適化問題を解くための準備として，東京証券取引所（東証）に上場している複数の金融資産（株式銘柄）の個別銘柄データを準備します．
今回はTOPIX 500に掲載されている銘柄から10個選び，投資対象の金融資産とします．これらの個別銘柄データはpandas-datareaderを使ってstooqから取得できます．データを取得したら，パネルデータを作りpandas.Dataframeオブジェクトとして整理します．

以下でパネルデータを作成するための関数を定義します．

def get_stockvalues(stockcode, start, end, use_ratio=False):
    # Get index data from https://stooq.com/
    df = stooq.StooqDailyReader(f"{stockcode}.jp", start, end).read()
    df = df.sort_values(by='Date',ascending=True)

    # Convert "Close" value on return ratio based on "Close" value
    if use_ratio:
        df = df.apply(lambda x: (x - x[0]) / x[0] )
    return df

def get_paneldata(stockcodes, start, end, use_ratio=False):
    # Use "Close" value only 
    dfs = []
    for sc in stockcodes:
        df = get_stockvalues(sc, start, end, use_ratio)[['Close']]
        df = df.rename(columns={'Close': sc})
        dfs.append(df)
    df_concat = pd.concat(dfs, axis=1)
    return df_concat

それでは実際に銘柄を選択し，パネルデータを作成します．

stockcodes=["2413", "2432", "3092", "3632", "3769", "4689", "4751", "4755", "6098", "9843"]

今回は上記10個の銘柄を選択しました．投資開始日を2019年1月1日と仮定し，2022年7月31日時点でのデータを取得します．

start = datetime.datetime(2019, 1, 1)
end = datetime.datetime(2022, 7, 31)
stockcodes=["4751", "2413", "2432", "4755", "3092", "3632", "9843", "3769", "4689", "6098"]

df = get_paneldata(stockcodes, start, end, use_ratio=True)
display(df) # Show dataframe
------
	4751	2413	2432	4755	3092	3632	9843	3769	4689	6098
Date										
2019-01-04	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
2019-01-07	0.034105	0.077451	0.025385	0.053719	0.086329	0.036361	0.014045	0.063335	0.029329	0.055017
2019-01-08	0.032885	0.064424	0.030355	0.078512	0.094271	0.043210	0.025280	0.084444	0.047640	0.063155
2019-01-09	0.058469	0.151473	0.043598	0.099174	0.170023	0.052299	0.060042	0.122222	0.069597	0.118176
2019-01-10	0.017053	0.143242	0.042497	0.114325	0.156772	0.050034	0.037570	0.162221	0.065951	0.116039
...	...	...	...	...	...	...	...	...	...	...
2022-07-25	0.345171	1.708894	0.110259	-0.100551	0.606797	1.035267	0.039530	1.360157	0.730807	0.853217
2022-07-26	0.342195	1.644725	0.104012	-0.121212	0.587425	1.040153	0.029828	1.380406	0.740319	0.847303
2022-07-27	0.346163	1.641965	0.100605	-0.125344	0.573034	0.986401	-0.004308	1.375906	0.741905	0.841388
2022-07-28	0.268786	2.009729	0.109691	-0.121212	0.573587	1.042596	-0.012572	1.400655	0.795410	0.921826
2022-07-29	0.310451	2.184296	0.114235	-0.097796	0.584657	1.062143	0.010425	1.470402	0.859221	0.952976
------

fig = plt.figure(figsize=(12,8))
plt.plot(df.index, df.values, label=df.columns)
plt.title("Trends in the Returns of Multiple Stocks", fontsize=18)
plt.xlabel("Date", fontsize=18)
plt.ylabel("Return", fontsize=18)
plt.legend(fontsize=18)
plt.tick_params(labelsize=18)
plt.show()
fig.savefig("data.pdf")

ここで，use_ratio=Trueとすることにより前処理として，各日付時点における終値を初日の終値ベースの収益率に変換しています．
以上によりデータフレームが作成され，下図のように各銘柄ごとの収益率の推移がプロットできました（凡例は各銘柄を示しています）．

Codable Model Optimizerを用いて最適化

それでは実際にPythonの最適化フレームワークである「Codable Model Optimizer」を用いて，このデータフレームで定義された二次計画問題（QP）を解き，ポートフォリオの最適な投資比率を求めていきます．

まず，対象銘柄のパネルデータdfから，必要な統計量を計算します．

ポートフォリオ内の銘柄間の共分散行列 $\Sigma$：

df.cov() # Covariance matrix
-----
        2413	    2432	    3092	    3632	    3769	    4689	    4751        4755	    6098	    9843
2413	2.751460	0.094242	0.498758	0.192367	1.032176	0.856012	0.604603	0.185212	0.506950	0.319627
2432	0.094242	0.022524	0.033668	0.019898	0.036389	0.019982	0.031072	0.019857	0.033207	0.004880
3092	0.498758	0.033668	0.169810	0.113825	0.230558	0.192363	0.156797	0.049032	0.175889	0.055580
3632	0.192367	0.019898	0.113825	0.214011	0.137703	0.126672	0.090277	-0.004010	0.156493	-0.001461
3769	1.032176	0.036389	0.230558	0.137703	0.458460	0.358120	0.262881	0.081857	0.253478	0.118675
4689	0.856012	0.019982	0.192363	0.126672	0.358120	0.365520	0.205512	0.053383	0.227503	0.113690
4751	0.604603	0.031072	0.156797	0.090277	0.262881	0.205512	0.181702	0.058047	0.172272	0.066235
4755	0.185212	0.019857	0.049032	-0.004010	0.081857	0.053383	0.058047	0.056955	0.052818	0.026389
6098	0.506950	0.033207	0.175889	0.156493	0.253478	0.227503	0.172272	0.052818	0.240297	0.054326
9843	0.319627	0.004880	0.055580	-0.001461	0.118675	0.113690	0.066235	0.026389	0.054326	0.052601

ポートフォリオ内の各銘柄の期待収益率 ${\bf r}$：

df.mean() # Expected returns
-----
2413    2.388957
2432    0.045658
3092    0.479059
3632    0.436813
3769    1.415249
4689    0.923626
4751    0.437981
4755    0.427737
6098    0.761135
9843    0.254608

ここから，「Codable Model Optimization」を使っていきます．

Problemオブジェクトの生成

はじめに，Problemオブジェクトを生成します．生成時に，引数を用いて，目的関数の最大化/最小化問題の設定を行います．is_max_problem=Trueとした場合は，最大化問題，is_max_problem=Falseとした場合は，最小化問題となります．今回は，リスクの最小化問題であるため，is_max_problem=Falseとします．

problem = Problem(is_max_problem=False)

変数の定義

最適化で利用する変数を定義します．「Codable Model Optimizer」は、整数型・連続値型・カテゴリ型を利用することができ，整数・連続値型では，上界・下界値を設定することができます．今回は投資比率を求めるため，連続値型を選択し，上界値「1」・下界値「0」に設定します．また，各銘柄ごとに投資比率を決定する必要があるため選択する銘柄の数だけ変数（$x$）を定義します．

x = [DoubleVariable(name=sc, lower=0, upper=1) for sc in stockcodes]
len(x) # The number of stoks.
-----
10

目的関数の設定

本問題はポートフォリオのリスクを最小化することが目的となりますので，目的関数をPythonの関数と引数のマッピング情報によって以下のように定義できます．今回は，10個の要素からなる1次元配列の変数と，2次元配列の共分散行列を引数として設定します．

# Define the objective function
def objective_function(var_x, para_cov):
    var_x = np.array(var_x)
    return np.dot(var_x.T, np.dot(para_cov, var_x))

problem += Objective(objective=objective_function,
                     args_map={'var_x': x, 'para_cov': cov})

また，引数のマッピング情報にVariableオブジェクトを設定した場合は，Variableオブジェクトの値が渡されます．Variableオブジェクトは，1~3次元配列にしてまとめて渡すことも可能です．

制約式の設定

制約式は，Variableオブジェクトを含めた不等式をProblemオブジェクトに加えることで定義できます．今回のポートフォリオ最適化問題では2つの制約式を追加します．1つ目は不等式制約としてポートフォリオの期待利益率（期待リターン）が目標利益率（目標リターン）以上となる制約，2つ目は投資比率の和が1となる制約を加えます．

tret = mean.mean()
problem += np.sum(mean * x) >= tret
problem += np.sum(x) == 1
print(tret)
------
0.7570822492119635

本来，目標リターンは自分で設定しますが今回は目標リターンとして，とりあえず期待リターンの平均値を用います．ポートフォリオ最適化問題においては金融資産を増やすことが目的となりますので，基本的にはこの目標リターンを1以上に設定することが望ましいでしょう．

最適化の実行

最適化を実行するには，Solverオブジェクトを生成，利用する最適化手法のMethodオブジェクトを引数としてsolve関数に設定する必要があります．最適化手法にはデフォルトとしてPenaltyAdjustmentMethodが提供されており，この手法は重み付き局所探索法(WLS; Weighting Local Search)という手法の一種で，与えられた制約のもとで，目的関数$f(x)$を最小化する変数$x$を探索する問題を考えます．詳細についてはこちらを参照してください．

solver = OptSolver(round_times=2)
method = PenaltyAdjustmentMethod(steps=40000)
answer, is_feasible = solver.solve(problem, method)
-----
answer:{'2413': 0.0, '2432': 0.016912671421771774, '3092': 0.0, '3632': 0.07716465883778946, '3769': 0.3669398519728481, '4689': 0.0031589825848488485, '4751': 0.004625689077339748, '4755': 0.3765422167661788, '6098': 0.001860173741888252, '9843': 0.15279575559733388}
answer_is_feasible:True

実行結果としては，“獲得された解"と"実行可能解フラグ"が返ってきます．獲得された解は辞書型で各変数の値が各のされています．keyは変数生成時に定義した各銘柄の'name',valueは各銘柄の投資比率となっています．最適解が得られていた場合，この解は期待リターンが自ら設定した目標リターン以上であるという条件のもとでリスクが最小となる投資比率を表しています．

効率的フロンティアを描く

上記の最適化の際には，目標リターンとして期待リターンの平均値を用いていましたが，本来は自分が求める目標リターンを定義する必要があります．ただし，当然のことながら目標リターンを大きくすればするほどリスクは大きくなります．このような関係をトレードオフと呼びます．

効率的フロンティアとは，このトレードオフの関係の上で最適なリスクとリターンの組み合わせの集合を指します．ここでは詳細な説明を省き，とりあえず効率的フロンティアを図示していきたいと思います．

目標リターンを複数パターン準備

trets = np.linspace(0, 1.5, 50)
print(trets)
-----
[0.         0.03061224 0.06122449 0.09183673 0.12244898 0.15306122
 0.18367347 0.21428571 0.24489796 0.2755102  0.30612245 0.33673469
 0.36734694 0.39795918 0.42857143 0.45918367 0.48979592 0.52040816
 0.55102041 0.58163265 0.6122449  0.64285714 0.67346939 0.70408163
 0.73469388 0.76530612 0.79591837 0.82653061 0.85714286 0.8877551
 0.91836735 0.94897959 0.97959184 1.01020408 1.04081633 1.07142857
 1.10204082 1.13265306 1.16326531 1.19387755 1.2244898  1.25510204
 1.28571429 1.31632653 1.34693878 1.37755102 1.40816327 1.43877551
 1.46938776 1.5       ]

まずは制約条件である目標リターンのパターンを網羅的に準備します．
今回は0から1.5までの50個を準備しました．

各目標リターンに対して最適化を実行

vols = []
rets = []
for i, tret in enumerate(trets):
    problem = Problem(is_max_problem=False)
    x = [DoubleVariable(name=sc, lower=0, upper=1) for sc in stockcodes]
    problem += Objective(objective=objective_function,
                     args_map={'var_x': x, 'para_cov': cov})
    problem += np.sum(x) == 1
    problem += np.sum(mean * x) >= tret
    solver = OptSolver(round_times=2)
    method = PenaltyAdjustmentMethod(steps=80000)
    answer, is_feasible = solver.solve(problem, method)
    res = list(answer.values())
    vols.append(np.sqrt(objective_function(res, cov)))
    rets.append(np.sum(mean * res))
vols = np.array(vols)
rets = np.array(rets)

続いて，各目標リターンのパターンごとに最適化を行い，最終的に獲得された解（投資比率）からリスクと期待リターンを算出し，その組み合わせを取得します．

リスク・リターンの関係性をプロット

得られた結果を図示していきます．

cm = plt.cm.get_cmap('RdYlBu')
fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(1, 1, 1)
mappable = ax.scatter(vols, rets, c=rets/vols, vmin=0, vmax=max(rets/vols), marker='o', s=50, cmap=cm)
plt.grid(True)
plt.tick_params(labelsize=18)
fig.colorbar(mappable, ax=ax)
plt.title("Efficient frontier", fontsize=18)
plt.xlabel("Expected volatility", fontsize=18)
plt.ylabel("Expected return", fontsize=18)
plt.show()

最適化が正確に実行されていれば，目標リターンのパターンに対して，それぞれ期待リターンが目標リターン以上である条件のもとでリスクが最小となる投資比率が求められます．まさに，これらの組み合わせの集合のことを効率的フロンティアと呼びます．この効率的フロンティアを描くことができれば，リスクとリターンのトレードオフ関係を把握することができ，より自分の要求にそった投資比率を決定することができます．
また，図のカラーバーは「期待リターン/リスク」の値を表しており，その値が大きくなるほどマーカーの色が青く，小さくなるほど赤く塗られています．ポートフォリオ最適化の文脈ではできるだけ少ないリスクで大きなリターンを得たいことが多いため，効率的フロンティアの上であっても，より青色で表される部分を選ぶことが理想と言えるでしょう．

補足〜効率的フロンティアが描けていない場合〜

目標リターンの設定が大きすぎる場合などには最適化フロンティアが正しく形成できないことがあります．
以下のように目標リターンを設定した場合の最適化結果をみてみましょう．

trets = np.linspace(0, 3.0, 50)
print(trets)
-----
[0.         0.06122449 0.12244898 0.18367347 0.24489796 0.30612245
 0.36734694 0.42857143 0.48979592 0.55102041 0.6122449  0.67346939
 0.73469388 0.79591837 0.85714286 0.91836735 0.97959184 1.04081633
 1.10204082 1.16326531 1.2244898  1.28571429 1.34693878 1.40816327
 1.46938776 1.53061224 1.59183673 1.65306122 1.71428571 1.7755102
 1.83673469 1.89795918 1.95918367 2.02040816 2.08163265 2.14285714
 2.20408163 2.26530612 2.32653061 2.3877551  2.44897959 2.51020408
 2.57142857 2.63265306 2.69387755 2.75510204 2.81632653 2.87755102
 2.93877551 3.        ]

上図の結果のように，目標リターンが大きすぎるとそもそも制約充足解が求められず，正しく効率的フロンティアが描けない場合もありますので注意しましょう．

また，目標リターンが十分大きすぎない場合であっても最適解が求められなければ効率的フロンティアは描くことができません．
そこで，以下のようにPenaltyAdjustmentMethodのstep数を少なくして最適化を実行してみましょう（目標リターンのパターンは0から1.5までの50個を準備）．

method = PenaltyAdjustmentMethod(steps=10000)

このように最適化手法の精度が悪いと最適解が求められないため，赤線部のように期待リターンに対して十分に小さいリスクを示す解が得られなくなってしまいます．このような場合についても正しく効率的フロンティアが描けないので注意しましょう．

最後に

今回，最適化フレームワーク「Codable Model Optimizer」でポートフォリオ最適化を行った結果を紹介しました．世の中には最適化ライブラリがたくさんありますが，本ライブラリは非常に扱いやすいライブラリだなと感じたので，ぜひ本記事を参考にして利用してみていただければ幸いです．また，今回参考にさせていただいた記事や，「Codable Model Optimizer」の紹介記事を改めて下記に記します．最後まで読んでいただきありがとうございました．

参考

免責事項

本内容はポートフォリオ理論や数理計画問題の一般論であり，投資を勧めるものでは一切ありません．
本投稿で計算した数値は仮定のものであり，その数値に一切責任は負いません．
本内容は個人の見解であり，所属組織とは一切関係ありません．