👋

ポートフォリオ最適化:マルコビッツモデルと進化計算(進化戦略)

2024/05/04に公開

はじめに

本記事ではポートフォリオ最適化の過程でマルコビッツモデルにより算出される投資のリスクを進化計算を用いて最適化することを目的とする.

金融に関しては素人のため,間違い,ご指摘等ございましたら,よろしくお願いいたします.

マルコビッツモデルについては下記記事を参考にした.

参考記事:ポートフォリオ最適化問題

以下で本記事の構成について説明する.

第1章では最適化とは何かについて具体的な例を交えながら説明し,本章の最後にマルコビッツモデルについて説明を行い,定式化を行う.

第2章では進化計算について今回用いるアルゴリズムの説明を行い,Pythonでの実装及び,引用元を示す.

第3章では実際にマルコビッツモデルにより算出される投資のリスクの最適化にチャレンジしてみる.

第4章では投資におけるリスクとリターンの関係性に着目する.

第5章で本記事のまとめを行う.

Appendixにはその他必要な情報を記載する.(データの収集コードなど)

第1章 最適化とは

本章では初めに最適化に関してざっくりとした説明を行う.
次に具体例を示す.
最後に本記事で扱うマルコビッツモデルの定式化と最適化の対象について説明を行う.

最適化のざっくりとした説明

皆さんは最適化という言葉を聞いたことがあるだろうか.最適化とはズバリ,何かを最小化(または最大化)するという行為そのものである.
日常生活において,皆さんが趣味に没頭する,ドライブに行く,ショッピングを楽しむ,筋トレをするなどといった行為も,自らの幸福を最適化(最大化)していると言えますね.
もちろん制約条件は自らの予算内です!!(借金は良くないです.)

なお,一般的に最適化では最小化を考えます.
maxmize x \leftrightarrow minimize -x (x\in\mathbb{R})
幸福を最大化 \leftrightarrow 不幸を最小化
などと考えればわかりやすいですね.

以下ではもう少し,実用的な具体例を示す.

具体例1 ポートフォリオ最適化

株式に投資するとき,一つの銘柄に全ての資源を投下するのではなく,
A社株を10%, B社株を15%,C社株を10%,...
のように投資の対象を分散させることで投資のリスクを低減させることが出来ます.

一般的にポートフォリオ最適化ではリスクとリターンのトレードオフを考慮しながら最適な資源の配分を決定します.

リスクを最小化する,ある程度のリスクを仮定したうえでリターンを最大化するなどといった最適化問題に帰着されます.

参考記事:ポートフォリオ最適化とは

具体例2 機械学習・AI

近年ブームになっているAIは最適化とは切っても切り離せない関係であることをご存知だろうか.
下記の図ニューラルネットワークが入力値x_m\in\mathbb{R}^n(下の図ではn = 2)に対して出力値\hat{y_m} \in \mathbb{R}を出力している様子である.
image.png

そしてニューラルネットワークを訓練するとき,入力値x_mに対する予測\hat{y_m}と実際の正解ラベルy_mとの損失関数(平均2乗誤差など)が小さくなるように各重みw_iを調整している.
損失関数の例:平均2乗誤差関数 MAE = \frac{1}{m} \Sigma_{i = 1} ^ m (\hat{y}_i - y_i)^2

すなわち,AIの訓練とはある意味,損失関数の最適化(\hat{y}yの誤差の最小化)であると考えられる.
(実際にはデータをミニバッチに分割し,損失関数による誤差を計算したのちに,誤差逆伝播法により各パラメータの勾配を求め,重みw_i,バイアスbなどを更新しています.)

マルコビッツモデル

マルコビッツモデルは変数\mathbf{x}=(銘柄1の購入割合,銘柄2の購入割合,..., 銘柄nの購入割合 ) \in \mathbb{R}^n (nは銘柄の数)に対して,投資のリスクや,リターンの期待値を算出することが出来るモデルである.

マルコビッツモデルの式などの詳細については下記記事をご参照ください.

参考記事:マルコビッツモデルの定式化

本記事では,投資のリスクを最小にするような投資配分を求めることに挑戦する.

※収益率とは各株が割安なのか割高なのかを表す指標であり,株価 / EPSで表される.
(EPSは1株あたりの純利益)

以下にマルコビッツモデルのpythonでの実装を示す.

def Markowitz(x, ave_r = ave_r, df = df):
    #ave_r = [銘柄1の収益率の平均,銘柄2の収益率の平均,銘柄3の収益率の平均, ...]
    #df→各銘柄の収益率が列ごとに格納されている.(今回は10期分)

    #空売り禁止の制約チェック.
    if np.any(x < 0.0):
        return float("inf")
    
    # 合計が1になるように正規化
    sum_x = np.sum(x)  
    x = np.round(x / sum_x, decimals=8)  
    
    # 平均リターン
    ave_rp = np.dot(ave_r, x)  
    
    # 平均リターンが負の場合も制約違反とみなす.
    if ave_rp < 0.0:  
        return float("inf")
    
    devs = np.dot(df.values, x) - ave_rp 
    
    # ポートフォリオのリスクを計算
    risk = np.sum(devs ** 2) / len(ave_r) 
    
    return risk  

第2章 進化計算とは

本章では,はじめに進化計算について説明を行う.
次に本記事で用いるアルゴリズムDXNESについて説明を行う.
最後にその実装と具体的な最適化の挙動をGIFのアニメーションで示す.

進化計算とは

進化計算とは生物が環境に合わせて進化する様子からインスピレーションを受けた最適化手法であり,最適化の過程で試行錯誤を繰り返すことで最適解,よりよい解の発見を目指す手法である.
以下にメリット,デメリットを示す.

メリット

  • 微分不可能な(勾配情報が使用できない)目的関数(タスク)に対しても
    最適化を行うことができる.
  • 機械学習などと異なり,(基本的に)事前の学習データが不要である.
    (訓練データと正解ラベルなど)
  • 関数景観が分からなくても,各タスクに適用できる.

デメリット

  • タスクに応じたサンプルサイズのチューニング,探索の終了条件の決定が難しい.
  • 探索開始時の正規分布のパラメータの設定がユーザパラメータである.
  • 学習率などはあくまでも開発者が時間と労力をかけてチューニングした推奨値であり,理論的に証明されたものではない.

なお本記事では,DXNESのサンプルサイズ,探索の終了条件は決め打ちですが,実際には
正規分布のパラメータ,サンプルサイズをタスクに応じて適切に設定する必要があることにご注意ください.

参考記事:進化計算とは
参考記事:進化計算のウィキペディア

DXNES

今回は進化計算の中でも,進化戦略の一種であるDXNESを使用する.
DXNESの探索の過程を説明し,実際に最適化している様子をGIFで紹介する.
はじめに探索の手順を示す.

1. 探索空間上に正規分布を初期化する.
2. 正規分布に従い解\mathbf{x} \in \mathbb{R}^nを生成(サンプリング)し,各解\mathbf{x} \in \mathbb{R}^nに対して評価値\in \mathbb{R}を計算する.
3. 評価値の小さい解\mathbf{x} \in \mathbb{R}^nが多く存在する方向に正規分布を移動させる.

1.の手順で正規分布を初期化したのちに2.と3.を繰り返すことで最適化を行う.

DXNESを詳しく知りたい方は下記記事をご参照ください.

参考記事:DXNESとは

次に,実際の最適化の挙動をRastrigin関数を用いて示す.

Rastrigin関数の式と関数景観 (n = 2)

  • Rastrigin関数の式
    f(\mathbf{x}) = 10 n + \Sigma_{i = 1}^n ((x_i - opt)^2 - 10\ \mathrm{cos}(2\pi \ (x_i - opt)))
    今回は最適解opt = (3.0, 3.0)とした

  • Rastrigin関数の景観
    ヒートマップ上では青い部分ほど評価値が良い値となっている.

Rastrigin関数を最適化する際の挙動

下のアニメーションでは白点が正規分布の中心,黒点がサンプリングされる解を表している.
アニメーションを見ると左下の(-3.0, -3.0)\in \mathbb{R}^2に初期化された正規分布の中心が,解\mathbf{x}\in \mathbb{R}^2をサンプリングしながらより良い評価値を持つ方向へ移動し,正規分布が最適解(3.0, 3.0) \in \mathbb{R}^2収束していることわかる.

DXNESの実装

ここではDXNESのPythonによる実装を示す.(ここは飛ばしてもらって結構です.)
実装,パラメータは下記の引用元に従った.

引用元

class solution:
    def __init__(self, dim):
        self.x = np.zeros(dim)
        self.z = np.zeros(dim)
        self.dim = dim
        self.eval = 0.0
            
    def __str__(self):
        sol = "個体:["
        for i in range(self.dim):
            if not i == self.dim - 1:
                sol += str (self.x[i]) + ", "
            else:
                sol += str (self.x[i]) + " ]\n"
                
        sol += f"評価値:{self.eval}"
        return sol
    
    def __lt__(self, other):
        return self.eval < other.eval
    
    def get_z(self):
        return self.z
    
    def set_z(self, z):
        self.z = z
    
    def get_x(self):
        return self.x

    def set_x(self, x):
        self.x = x
        
    def set_eval(self, eval):
        self.eval = eval
        
    def get_eval(self):
        return self.eval
    
class dxnes:
    def __init__(self, dim, pop_size, m, sigma, seed):
        np.random.seed(seed)
        #集団の保持
        self.population = []
        for _ in range(pop_size):
            sol = solution(dim)
            self.population.append(sol)
        #集団サイズ
        self.pop_size = pop_size
        #次元数
        self.dim = dim
        #平均ベクトル,ステップサイズ,正規化変換行列
        self.m = np.ones(dim) * m
        self.sigma = sigma
        self.B = np.eye(dim)
        
        #重みの初期化
        self.w = np.zeros(self.pop_size)
        
        #rankベースの重みの初期化
        self.w_hat = np.zeros(pop_size)
        self.w_rank = np.zeros(pop_size)
        for i in range(pop_size):
            tmp_value = np.log((float(pop_size) / 2.0) + 1.0) - np.log(i + 1)
            if 0.0 < tmp_value:
                self.w_hat[i] = tmp_value
                
        sum = float (self.w_hat.sum())        
        for i in range(pop_size):
            self.w_rank[i] = (self.w_hat[i] / sum) - (1.0 / float(pop_size))
            
        #distベースの重みの初期化
        self.w_dist = np.zeros(pop_size)
            
        #h^(-1)(n)の計算をニュートン法を用いて行う.
        a = 1.0
        for i in range(1000):
            tmp_1 = (1.0 + float(a ** 2)) * np.exp(float(a ** 2) / 2.0)
            tmp_1 = (tmp_1 / 0.24) - 10.0
            tmp_1 = tmp_1 - float(self.dim)
            tmp_2 = float (a) * float(a ** 2) * np.exp(float(a ** 2) / 2.0)
            tmp_2 = tmp_2 / 0.12
            a = a - float(tmp_1) / float(tmp_2)
        
        self.n_inv = a
        
        #alphaの計算を行う.
        self.alpha = 1.0 * self.n_inv
        if np.sqrt(float(pop_size)/ float(dim)) < 1.0:
            self.alpha = self.n_inv * np.sqrt(float(pop_size)/ float(dim))
            
        #学習率の設定を行う.
        #平均ベクトルの学習率
        self.Eta_m = 1.0
        
        #正規化変換行列の学習率
        self.Eta_B = 0.0
        self.Eta_B_move = (float(pop_size) + 2.0 * float(dim)) / (float(pop_size) + 2.0 * float(dim) * float(dim) + 100.0)
        
        if np.sqrt(float(pop_size)/ float(dim)) < 1.0:
            self.Eta_B_move = self.Eta_B_move * np.sqrt(float(pop_size)/ float(dim)) 
        
        self.Eta_B_stag = (float(pop_size)) / (float(pop_size) + 2.0 * float(dim) * float(dim) + 100.0)
        self.Eta_B_conv = (float(pop_size)) / (float(pop_size) + 2.0 * float(dim) * float(dim) + 100.0)
        
        #ステップサイズの学習率
        self.Eta_sigma = 0.0
        self.Eta_sigma_move = 1.0
        self.Eta_sigma_stag = 0.5 * (1.0 + (float(pop_size) / (float(pop_size) + 2.0 * float(dim))))
        self.Eta_sigma_conv = 2.0 * self.Eta_sigma_stag
        
        #c_sigmaの初期化
        mu_eff = 0.0
        for i in range(pop_size):
            mu_eff += (float(self.w_rank[i]) + (1.0 / float(self.pop_size))) * (float(self.w_rank[i]) + (1.0 / float(self.pop_size)))
        
        self.mu_eff = 1.0 / mu_eff
        
        tmp_1 = 1.0 / (2.0 * np.log(float(dim) + 1.0))
        tmp_2 = (self.mu_eff + 2.0) / (float(dim) + self.mu_eff + 5.0)   
        
        self.c_sigma = tmp_1 * tmp_2
        
        #進化パスの初期化
        self.p_sigma = np.zeros(dim)
        
        #ε(フェーズの判定に用いる閾値)→重み,学習率の決定に用いる
        self.epsilon = np.sqrt(dim) * (1.0 - (1.0 / (4.0 * float(dim))) + (1.0 / (21.0 * float (dim) * float (dim))))
        
        #自然勾配の計算に用いるパラメータ
        self.G_B = np.eye(dim)
        self.G_M = np.eye(dim)
        self.G_delta = np.zeros(dim)
        self.G_sigma = 0.0

        #計算用の行列,ベクトル
        self.vecter = np.zeros(dim)
        self.matrix = np.eye(dim)
        
    #個体のサンプリングを行う.     
    def sampling(self):
        count = int (self.pop_size / 2)
        for i in range(count):
            z = np.random.normal(0.0, 1.0, self.dim)
            self.population[2 * i].z = z
            self.population[2 * i].x = self.sigma *  np.real(np.dot(self.B, z)) + self.m
            
            z = -1.0 * z
            self.population[2 * i + 1].z = z
            self.population[2 * i + 1].x = self.sigma *  np.dot(self.B, z) + self.m
            #print(self.population[2 * i + 1].x)
            
    #集団内の個体をソートする
    def sort(self):
        self.population.sort()
            
    #進化パスの計算を行う.
    def calc_evolution_path(self):
        tmp_value = 0.0
        for i in range(self.pop_size):
            tmp_value += self.w_rank[i] * self.population[i].get_z()
        tmp_value *= np.sqrt(self.c_sigma * (2.0 - self.c_sigma) * self.mu_eff)
        tmp_value += (1.0 - self.c_sigma) * self.p_sigma
        self.p_sigma = tmp_value
        
    #distベースの重みw_distの計算
    def calc_dist_weight(self):
        tmp_value = 0.0
        for i in range(self.pop_size):
            tmp_value += self.w_hat[i] * np.exp(self.alpha * np.linalg.norm(self.population[i].get_z()))
        
        for i in range(self.pop_size):
            self.w_dist[i] = (((self.w_hat[i]) * np.exp(self.alpha * np.linalg.norm(self.population[i].get_z()))) / (tmp_value)) - (1.0 / float(self.pop_size))
        
        self.w = self.w_dist                                     
        
        
    #進化パスのノルムに応じた重みの設定を行う.
    def deside_weight(self):
        norm = np.linalg.norm(self.p_sigma)
        if norm >= self.epsilon:
            self.calc_dist_weight()
        else:
            self.w = self.w_rank
            
    #ステップサイズ,正規化変換行列の学習率を決める
    def decide_learning_rate(self):
        norm = np.linalg.norm(self.p_sigma)
        #move, #stag, #convの学習率の設定
        if norm >= self.epsilon:
            self.Eta_sigma = self.Eta_sigma_move
            self.Eta_B = self.Eta_B_move
        elif norm >= 0.1 * self.epsilon:
            self.Eta_sigma = self.Eta_sigma_stag
            self.Eta_B = self.Eta_B_stag
        else:
            self.Eta_sigma = self.Eta_sigma_conv
            self.Eta_B = self.Eta_B_conv
            
    #自然勾配の計算を行う.
    def calc_Natural_Gradient(self):
        #G_delta, G_Mの計算
        pop_size = len(self.population)
        eye_matrix = np.eye(self.dim)
        self.G_delta = np.zeros(self.dim)
        self.G_M = np.eye(self.dim) * 0.0
        for i in range(pop_size):
            self.G_delta += self.w[i] * self.population[i].get_z()
            self.G_M += self.w[i] * (np.dot(self.population[i].get_z().reshape(self.dim, 1), self.population[i].get_z().reshape(1, self.dim)) - eye_matrix)
        #G_sigma, G_Bの計算を行う.
        self.G_sigma = float(np.trace(self.G_M)) / float(self.dim)
        self.G_B = self.G_M - self.G_sigma * eye_matrix
        
    def update_parameters(self):
        #平均ベクトルmの更新
        self.vecter =  self.Eta_m * self.sigma * np.dot(self.B, self.G_delta)
        self.m = self.m + self.vecter
        
        #ステップサイズσの更新
        self.sigma = self.sigma * np.exp(self.Eta_sigma * self.G_sigma / 2.0)
        
        #正規化変換行列Bの更新
        self.matrix = self.Eta_B * self.G_B / 2.0
        eig, P = np.linalg.eig(self.matrix)
        self.matrix = np.linalg.inv(P) @ self.matrix @ P
        
        for i in range(self.dim):
            self.matrix[i][i] = np.exp(self.matrix[i][i])
        
        self.matrix = np.linalg.inv(P) @ self.matrix @ P
        
        self.B = np.dot(self.B, self.matrix)
        
    def do_one_generation(self):
        self.sort()
        self.calc_evolution_path()
        self.deside_weight()
        self.decide_learning_rate()
        self.calc_Natural_Gradient()
        self.update_parameters()
        
    def get_population(self):
        return self.population
    
    def get_best_solution(self):
        self.population.sort()
        return self.population[0]
    
    def get_worst_solution(self):
        self.population.sort()
        return self.population[self.pop_size - 1]
    
    def get_m(self):
        return self.m

第3章 マルコビッツモデルにより算出されるリスクの最適化

本章ではマルコビッツモデルにより算出されるリスクを最小化する資源配分を求める.
はじめに,用いるデータの説明を行う.
なお,株のデータの収集に用いたコードについては本記事の最後にAppendixとして記載する.
次に,最適化のコードを示す.
最後に,最適化の結果について考えてみる.

最適化ではDXNESを用いて以下の手順で行う.

  1. 正規分布を初期化する.
  2. 正規分布に従い,解\mathbf{x} \in \mathbb{R}^n(nは銘柄数)を生成(サンプリング)し,マルコビッツモデルによりリスクを計算する.
  3. リスクの小さい解\mathbf{x} \in \mathbb{R}^nが多く存在する方向へ正規分布のパラメータを更新する.

2.と3.を打ち切り世代数まで繰り返す.

用いるデータ

データ収集に用いたコードはAppendixに記載した.

企業名

今回はS&P500より下記の24社の企業を用いる.
24社の理由は筆者都合です.EPSの取得に用いたサービスalphavantageは無料版だと1日25回までしか利用できない.
当初は筆者の興味本位でS&P500ではないBTC-USD(ビットコイン)を含めた25社を用いようとしたが,BTC-USDだけEPSが取得できずに24社となった.

    'AAPL', 'MSFT', 'GOOGL', 'AMZN', 'TSLA', 'NVDA', 'INTC', 'NFLX', 'AMD',  
    'PYPL', 'V', 'DIS', 'CRM', 'JPM', 
      'MCD', 'ACN', 'SBUX',
    "ADBE", "CSCO", "IBM", "GS", "NKE", "RTX",
    "KO"

企業は下記サイトを参考にいたしました.
参考記事:S&P500

EPS情報

EPS情報は各銘柄に対して四半期ごとの決済情報を参考にした.下記日付を四半期の終了日とする10期分のEPS情報を用いる.

'2023-12-31', '2023-09-30', '2023-06-30', '2023-03-31', '2022-12-31', '2022-09-30', '2022-06-30', '2022-03-31', '2021-12-31', '2021-09-30'

株価情報

株価は原則として各四半期決済の終了日の終値を用いた.しかし,四半期決済の終了日の株価が取得できなかった部分については,終了日から最も近く,値が存在する日付の終値を用いた.具体的な日付は下記である.

'2023-12-29', '2023-09-29', '2023-06-30', '2023-03-31', '2022-12-30', '2022-09-30', '2022-06-30', '2022-03-31', '2021-12-31', '2021-09-30'

最適化のコード

次元数dimは今回対象にしている企業の数である.

def evaluate(func, poplation):
    pop_size = len(poplation)
    for i in range(pop_size):
        eval = func(poplation[i].get_x())
        poplation[i].set_eval(eval)

dim = 24
nes = dxnes(dim = dim, pop_size = 80, m = 0.5, sigma = 0.25, seed = 1)  
for itr in range(1000):
    nes.sampling()
    evaluate(func = Markowitz, poplation = nes.get_population())
    nes.do_one_generation()
    print(nes.get_best_solution())

#得られた最良個体を正規化し,100を掛けることで%表示に直す.
result = (nes.get_best_solution().get_x() / np.sum(nes.get_best_solution().get_x())) * 100.0
print(result)

最適化の結果

本節では最適化の結果について考えてみる.

なお,以下で示す投資配分は最適化の過程で得られた解\mathbf{x} \in \mathbb{R}^nの一つに過ぎない.
すなわち,これが真の最適解かどうかは誰にも分からない.

現在は,正規分布のパラメータ,サンプルサイズは決め打ちであるが,ここをチューニングすることでより良い解が得られるかもしれない.

最適解は別の配分かもしれないし,複数個あるかもしれない.

具体的に最適化の過程において,投資のリスクの推移を表示し,得られた投資配分に関する投資配分を示してみる.

最適化の過程におけるリスクの推移と投資配分

  • 最適化の過程におけるリスクの推移
    下記のグラフは
    1枚目のグラフが横軸が世代数,縦軸がその世代で最もリスクの低かった投資のリスクを表している.
    2枚目のグラフは横軸が,世代数,縦軸がその世代時点における最小のリスクを表している.
    1枚目のグラフの途切れている部分は,最良個体のリスクがinf,すなわち,その世代の集団内に実行可能な個体が存在しなかった世代である.
    2枚のグラフより世代数が進み,最終的に,リスクの小さな投資配分が見つかっていることが分かる.

  • 最適化された投資配分(の一例)
    下記の投資配分は1000世代目に得られた最良個体をもとに投資配分を算出した結果である.

第4章 投資におけるリスクとリターンの関係性

本章では投資におけるリスクとリターンの関係について考える.
はじめにリスクとリターンの2目的最適化について定式化を行う.
そして,リスクとリターンの関係性に着目する.
もちろん,投資においてリスクとリターンはトレードオフの関係である.
すなわち,リスクを小さくすれば,リターンも小さくなり,逆に大きなリターンを狙うとそれだけリスクも大きくなるためこの関係性を確認する.

リスクとリターンの2目的最適化

本節では2目的最適化について定式化を行う.
具体的に最適化の対象を下記のように修正する.
なお,maximize (リターン) \leftrightarrow minimize (-1.0 \times リターン) に注意.

  • 修正前
    リスク(第3章まで)
  • 修正後
    リスク \times w + (1.0 - w) \times リターン \times (-1.0)

ここで0.0\leq w \leq1.0であり,リスクの考慮度合いを表す.
例えば,

  • w = 1.0のとき(最大限リスクを考慮)
    リスク \times w + (1.0 - w) \times リターン \times (-1.0)
    = リスク
    となり,最適化の対象はリスクとなる.

  • w = 0.0のとき(リスクを無視している)
    リスク \times w + (1.0 - w) \times リターン \times (-1.0)
    = リターン \times (-1.0)
    となり,最適化の対象はリターンとなる.

  • w = 0.5のとき
    リスク \times w + (1.0 - w) \times リターン \times (-1.0)
    = リスク \times 0.5 + 0.5 \times リターン \times (-1.0)
    = 0.5(リスク + リターン\times (-1.0))

すなわち,リスクとリターンのどちらも同程度考慮していることになる.

ここで修正後のマルコビッツモデルの実装を示す.

def Markowitz(x, weight,  ave_r = ave_r, df = df):
    #ave_r = [銘柄1の収益率の平均,銘柄2の収益率の平均,銘柄3の収益率の平均, ...]
    #df→各銘柄の収益率が列ごとに格納されている.(今回は10期分)
    
    #空売り禁止の制約チェック.
    if np.any(x < 0.0):
        return float("inf"), float("inf"), float("inf")
    
    # 合計が1になるように正規化
    sum_x = np.sum(x)  
    x = np.round(x / sum_x, decimals=8)  
    
    # 平均リターン
    ave_rp = np.dot(ave_r, x)  
    
    # 平均リターンが負の場合も制約違反として大きな値を返す
    if ave_rp < 0.0:  
        return float("inf"), float("inf"), float("inf") 
        
    devs = np.dot(df.values, x) - ave_rp 
    
    # ポートフォリオのリスクを計算
    risk = np.sum(devs ** 2) / len(ave_r) 
   
    #重みを付けた評価値を計算
    eval = (ave_rp * (-1.0)) * (1.0 - weight) + risk * weight
    
    return eval, risk, ave_rp  

最適化の結果

本節では0.0\leq w \leq1.0(リスクの考慮度合い)を少しずつ変化させたときのリスクとリターンの関係を可視化してみる.以下にその結果を示す.

  • 横軸:w, 縦軸:リターン

  • 横軸:リスク,縦軸:リターン

上記の結果よりリスクを考慮するとリターンは小さくなり,リスクを考慮しないとリターンが大きくなることから,リスクとリターンがトレードオフの関係であることが分かった.

第5章 まとめ

本記事では,下記の流れで話を進めた.
-第1章 最適化とマルコビッツモデルの定式化
-第2章 進化計算の説明,挙動のアニメーション,実装(Python)
-第3章 マルコビッツモデルにより算出される投資のリスクの最適化.
-第4章 リスクとリターンの2目的最適化

なお本記事では,DXNESのサンプルサイズ,探索の終了条件は決め打ちであったが,実用的に使用する際は特に
サンプルサイズをタスクに応じて適切にチューニングすることにご注意ください.

筆者は最近,NEATに興味を持っているのでまた気が向いたら何かを書きたいと思います.

最後までお読みいただき,ありがとうございました.

感想,ご指摘などがあれば,ぜひよろしくお願いいたします.

Appendix

ここでは本記事で用いたEPS, 株価の取得コードを記載いたします.

EPSの取得に用いたコード

(BTC-USDが取得できずに,24社となった.)

import requests
import pandas as pd
def get_financial_data_EPS(symbols, api_key):
    dfs = []
    for symbol in symbols:
        # Alpha Vantage APIのエンドポイント
        endpoint = f'https://www.alphavantage.co/query?function=EARNINGS&symbol={symbol}&apikey={api_key}'

        # APIからデータを取得
        response = requests.get(endpoint)
        data = response.json()
        
        # 'quarterlyEarnings' キーが存在するか確認
        if 'quarterlyEarnings' in data:
            # DataFrameに変換
            df_EPS = pd.DataFrame(data['quarterlyEarnings'])
            
            # 企業名をカラム名に追加
            df_EPS.columns = [f'{col}_{symbol}' for col in df_EPS.columns]

            dfs.append(df_EPS.head(10))  # 最初の10期分のデータを追加

    # 複数の企業のデータを結合して返す
    df_combined = pd.concat(dfs, axis=1)
    return df_combined

# 企業シンボルとAPIキーを設定
us_symbols = [
    'AAPL', 'MSFT', 'GOOGL', 'AMZN', 'TSLA', 'NVDA', 'INTC', 'NFLX', 'AMD',  
    'PYPL', 'V', 'DIS', 'CRM', 'JPM', 
     'BTC-USD', 'MCD', 'ACN', 'SBUX',
    "ADBE", "CSCO", "IBM", "GS", "NKE", "RTX",
    "KO"                
]
api_key = '自分のAPIキー'

# EPSデータを取得
df_EPS_combined = get_financial_data_EPS(us_symbols, api_key)

#EPS情報を取り出す.
reported_EPS_cols = [col for col in df_EPS_combined.columns if 'reportedEPS' in col]
df_reported_EPS = df_EPS_combined[reported_EPS_cols]

#欠値を埋める操作を行っている.
import numpy as np
for i in range(10):
    for j in range(24):
        if df_reported_EPS.iloc[i, j] == "None":
            df_reported_EPS.iloc[i, j] = np.nan
#欠値に直後の値を埋める.
df_reported_EPS = df_reported_EPS.fillna(method = "bfill")

株価の取得に用いたコード

import yfinance as yf
import pandas as pd

def get_stock_prices(symbols, target_dates):
    dfs = []
    for symbol in symbols:
        stock_data = yf.download(symbol, start='2021-01-01', end='2024-02-02')
        end_prices = stock_data['Close']
        end_prices.name = f'End Price_{symbol}'
        dfs.append(end_prices)

    df_combined = pd.concat(dfs, axis=1)
    df_filtered = df_combined.loc[target_dates]
    return df_filtered

us_symbols = [
    'AAPL', 'MSFT', 'GOOGL', 'AMZN', 'TSLA', 'NVDA', 'INTC', 'NFLX', 'AMD',  
    'PYPL', 'V', 'DIS', 'CRM', 'JPM', 
      'MCD', 'ACN', 'SBUX',
    "ADBE", "CSCO", "IBM", "GS", "NKE", "RTX",
    "KO"]


target_dates = ['2023-12-29', '2023-09-29', '2023-06-30', '2023-03-31', '2022-12-30', '2022-09-30', '2022-06-30', '2022-03-31', '2021-12-31', '2021-09-30']

df_stock_prices = get_stock_prices(us_symbols, target_dates)
df_stock_prices

Discussion