マルコビッツモデルと差分進化法(Differential Evolution)

はじめに

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

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

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

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

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

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

第2章では進化計算について今回用いるアルゴリズムである差分進化法(DE)の説明を行い，疑似コード，Pythonでの実装及び，参考記事を示す．

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

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

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

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

第1章 最適化とは

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

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

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

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

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

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

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

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

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

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

マルコビッツモデル

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

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

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

本記事では，マルコビッツモデルにより算出される投資のリスクを最小にするような投資配分
$\mathbf{x}$ $\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期分なので10×銘柄数のDataFrame)

    #空売り禁止の制約チェック．
    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章 進化計算とは

本章では，はじめに進化計算について説明を行う．
次に本記事で用いるアルゴリズムである差分進化法(DE)について説明を行う．
次に，DEの疑似コードとPythonの実装を示す．
最後にベンチマーク関数を用いて，DEの具体的な最適化の挙動をGIFのアニメーションで示す．

進化計算とは

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

メリット

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

デメリット

• タスクに応じたサンプルサイズ（集団内の個体数）のチューニング，探索の終了条件の決定が難しい．
• ユーザパラメータがそこそこ多く，探索の性能にも割と影響する．

参考記事：進化計算のウィキペディア

差分進化法(Differential Evolution)

今回は進化計算の中でも，差分進化法（以下DEと呼ぶ）を使用する．
DEは多くの分野で応用されている最適化手法の一種である．
DEの探索の過程を説明し，実際に最適化している様子をGIFで紹介する．
はじめに探索の手順を示す．
1.集団を探索空間上に初期化し，集団内の各個体に評価値を与える．
2.集団内の個体を用いて子個体を生成し，子個体に評価値を与える．
3.親個体（集団内の個体）と子個体を比較し，より最適解に近い（より評価値の小さい）個体を次世代の集団に残す．
4. 2と3を探索終了時まで繰り返す．

具体的な探索過程は後述のベンチマーク関数の一種であるRastrigin関数における最適化のGIFのアニメーションを参照してイメージを付けていただきたい．

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

引用元：DEの挙動，歴史

引用元：DEのウィキペディア

参考記事：ベンチマーク関数

次に，実際の最適化の挙動を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)$に近い解が次世代に引き継がれることを繰り返し，最適解$(3.0, 3.0)$を発見していることが分かる．

DEの実装

ここでは筆者が作成したDEの疑似コードとPythonによる実装を示す．

DEの疑似コード

DEのPythonによる実装

class solution:
    def __init__(self, dim):
        self.vector = np.zeros(dim)
        self.aim = np.zeros(dim)
        self.eval = 0.0
        self.eval_aim = 0.0
        self.dim = dim
            
    def __str__(self):
        sol = "個体：["
        for i in range(self.dim):
            if not i == self.dim - 1:
                sol += str (self.vector[i]) + ", "
            else:
                sol += str (self.vector[i]) + " ]\n"
                
        sol += f"評価値：{self.eval}"
        return sol
    
    def __lt__(self, other):
        return self.eval < other.eval
    
    def get_vector(self):
        return self.vector
    
    def set_vector(self, x):
        self.vector = x
        
    def get_aim(self):
        return self.aim
        
    def set_aim(self, x):
        self.aim = x
        
    def get_eval(self):
        return self.eval
    
    def set_eval(self, eval):
        self.eval = eval
    
    def get_eval_aim(self):
        return self.eval_aim

    def set_eval_aim(self, eval):
        self.eval_aim = eval
    
class DE:
    #スケーリングパラメータ
    F = 0.5
    #交叉率
    CR = 0.75
    
    def __init__(self, pop_size, dim, seed, max, min):
        #シード値の設定
        np.random.seed(seed)
        
        #集団サイズ
        self.pop_size = pop_size
        
        #最大値と最小値
        self.max = max
        self.min = min
        
        #交叉に用いるインデックスを求めるのに使うリスト
        self.all_list = []
        for i in range(self.pop_size):
            self.all_list.append(i)
        
        #集団の保持
        self.population = []
        for _ in range(pop_size):
            sol = solution(dim)
            self.population.append(sol)
        
        #次元数
        self.dim = dim
        
    def initialize_population(self, function):
        for i in range(self.pop_size):
            #ベクトルの生成と評価値の計算
            vec = np.random.rand(self.dim) * (self.max - self.min) + self.min
            eval = function(vec)
            
            #ベクトルと評価値をそれぞれ個体にセットする．
            self.population[i].set_vector(vec)
            self.population[i].set_eval(eval)
            
    def generate_aim_vector(self, function):
        
        for i in range(self.pop_size):
            cross_index = np.random.randint(0, self.dim)
            
            select_list = copy.deepcopy(self.all_list)
            select_list.pop(i)
            
            #交叉に用いる親子体の選択
            index_list = np.random.choice(select_list, 3, replace = False)
            parent_0 = self.population[index_list[0]].get_vector()
            parent_1 = self.population[index_list[1]].get_vector()
            parent_2 = self.population[index_list[2]].get_vector()
            
            #交叉を施したベクトルの作成
            mutaion_vector = parent_0 + DE.F * (parent_1 - parent_2)
            
            #比較対象のベクトルを作成する．
            aim_vec = np.zeros(self.dim)
            for j in range(self.dim):
                if np.random.rand() < DE.CR or j == cross_index:
                    aim_vec[j] = mutaion_vector[j]
                else:
                    aim_vec[j] = self.population[i].get_vector()[j]
            
            evel_aim = function(aim_vec)
        
            #aim_vectorとaim_evalをセットする．
            self.population[i].set_aim(aim_vec)
            self.population[i].set_eval_aim(evel_aim)
    
    def update_population(self):
        #aimと解の比較を行い，集団を更新する．
        for i in range(self.pop_size):
            #xとaimの評価値同士を比較し，aimの方が良ければpopulationの個体として更新．
            if self.population[i].get_eval_aim() < self.population[i].get_eval():
                self.population[i].set_vector(self.population[i].get_aim())
                self.population[i].set_eval(self.population[i].get_eval_aim())
                
    def do_one_generation(self, function):
            self.generate_aim_vector(function)
            self.update_population()
                
    def get_best_solution(self):
        return sorted(self.population)[0]
    
    def get_sorted_solutions(self):
        return sorted(self.population)

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

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

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

1. 集団を初期化する．
2. 集団より子個体を生成し，マルコビッツモデルによりリスクを計算する．
3. 親個体（集団内の個体）と子個体の評価値（リスク）を比較し，集団を更新する．

2.と3.を探索終了まで繰り返す．

用いるデータ

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

企業名

今回はS＆P500より下記の24社の企業を用いる．
24社の理由は筆者都合です．EPSの取得に用いたサービスalphavantageは無料版だと1日25回までしか利用できない．

alphavantage

当初は筆者の興味本位で25社のうち1社だけ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は今回対象にしている企業の数である．

de = DE(pop_size = 200, dim = 24, max = 1.0, min = 0.0, seed = 0)

function = Markowitz
de.initialize_population(function = function)
for i in range(1000):
    de.do_one_generation(function = function)
    print(de.get_best_solution())
    
#最良個体を正規化し，100%表示に直す
result_1 = de.get_sorted_solutions()[0].get_vector() / np.sum(de.get_sorted_solutions()[0].get_vector()) * 100.0

最適化の結果

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

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

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

具体的に最適化の過程において，投資のリスクの推移を表示し，得られた解$\mathbf{x} \in \mathbb{R}^n$を用いて投資配分を示してみる．

具体的に集団内200個体の中の最良個体(1st)，2番目の個体(2nd)，3番目の個体(3rd)，中間順位(mid)の個体のリスクの推移を示す．

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

• 最適化の過程におけるリスクの推移
  横軸が世代数，縦軸がリスクを表している．(表示は250世代までにしている．)
  

• 最適化された投資配分（の一例）
  下記の投資配分は各々，5,25,50,100,150,200,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目的最適化

気が向いたらまた何か書こうと思います．

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

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

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