QuantXでBlack-Litterman Modelを実装してみた

# はじめに
[株式会社Smart Trade](https://smarttrade.co.jp/)が提供している投資アルゴリズム開発プラットフォーム[QuantX Factory](https://quantx.io/)で，ポートフォリオ選択についての数理モデルであるBlack-Litterman Model[1][2][3][4]を実装しました．

# Black-Litterman Modelとは
現代ポートフォリオ理論は，1952年にHarry Markowitzが投資収益率の分布についての平均と分散のみに着目した平均分散分析（mean-variance analysis）に関する論文[5]を発表してからその歴史が始まります．

その後，1964年にWilliam F. Sharpeによって現代ポートフォリオ理論から発展したモデルである資本資産価格モデル（capital asset pricing model），通称CAPM[6]が発表されました．

この記事で扱うBlack-Litterman Modelは1992年にFischer BlackとRobert Littermanが上記の2つの現代ポートフォリオ理論の中核に位置する理論を実践するに当たって出くわす問題を克服したモデルです．

## その特徴
Markowitzの平均分散法には，一部の資産に偏った配分となる，推定に用いるデータ期間や目標リターンを少し変更しただけで極端に異なる配分結果となるなど，実務上の問題があります．[^1]

また，資産の共分散は適切に推定されているものの，期待リターンのもっともらしい推定値を導くのは難しいという問題があります．

この問題を期待リターンの推定値を必要としないことで解決しました．そのかわり，当初の期待リターンは均衡における資産配分が市場で観測されるものと同じとなるような期待リターンであると仮定します．

期待リターンが市場で観測されるリターンとどれほど違うのかと，代替的な仮定をどれほど信用するのかの程度のみが必要となります．このため，望ましい（平均分散的に効率的な）資産配分が計算可能になります．

一般に，ポートフォリオについての制約がある時，例えば，空売りが許容されない時，最適なポートフォリオを組むもっとも簡単な方法は資産の期待リターンを作るためにBlack-Litterman Modelを用いて，平均分散分析によって制約つき最適化問題を解くことです．[^2]

# 実装

<details>
  <summary>コードの全体像</summary>
  <div>

```python
import pandas as pd
import numpy as np

def initialize(ctx):
  ctx.configure(
    channels={
      "jp.stock": {
        "symbols": [
          "jp.stock.2914", #日本たばこ産業
          "jp.stock.3382", #セブン&アイ・ホールディングス
          "jp.stock.4063", #信越化学工業
          "jp.stock.4452", #花王
          "jp.stock.4502", #武田薬品工業
          "jp.stock.4503", #アステラス製薬
          "jp.stock.6098", #リクルートホールディングス
          "jp.stock.6501", #日立製作所
          "jp.stock.6752", #パナソニック
          "jp.stock.6758", #ソニー
          "jp.stock.6861", #キーエンス
          "jp.stock.6954", #ファナック
          "jp.stock.6981", #村田製作所
          "jp.stock.7203", #トヨタ自動車
          "jp.stock.7267", #本田技研工業
          "jp.stock.7751", #キヤノン
          "jp.stock.7974", #任天堂
          "jp.stock.8031", #三井物産
          "jp.stock.8058", #三菱商事
          "jp.stock.8306", #三菱UFJフィナンシャル・グループ
          "jp.stock.8316", #三井住友フィナンシャルグループ
          "jp.stock.8411", #みずほフィナンシャルグループ
          "jp.stock.8766", #東京海上ホールディングス
          "jp.stock.8802", #三菱地所
          "jp.stock.9020", #東日本旅客鉄道
          "jp.stock.9022", #東海旅客鉄道
          "jp.stock.9432", #日本電信電話
          "jp.stock.9433", #KDDI
          "jp.stock.9437", #NTTドコモ
          "jp.stock.9984", #ソフトバンクグループ
        ],
        "columns": [
          #"open_price_adj", # 始値(株式分割調整後)
          #"high_price_adj", # 高値(株式分割調整後)
          #"low_price_adj",  # 安値(株式分割調整後)
          #"volume_adj",     # 出来高
          #"txn_volume",     # 売買代金
          "close_price_adj", # 終値(株式分割調整後)
          "market_capitalization", # 時価総額(概算)
        ]
      }
    }
  )
  
  def _my_signal(data):
    cp = data["close_price_adj"].fillna(method="ffill")
    mkt_cap = data["market_capitalization"].fillna(method="ffill")
    
    sum_mkt_cap = pd.Series(data=0, index=cp.index)
    weight = pd.DataFrame(data=0, index=cp.index, columns=[])
    daily = pd.DataFrame(data=0, index=cp.index, columns=cp.columns)
    df_w_posterior = pd.DataFrame(data=0, index=cp.index[245:], columns=cp.columns)
    
    for sym, val in cp.items():
      sum_mkt_cap += mkt_cap[sym]
    
    for sym, val in cp.items():
      weight[sym] = mkt_cap[sym] / sum_mkt_cap
      
    # 2x30 (投資家のビュー)
    P = np.array(
          [
            [-0.5,-0.25,0,0,0,0,0,-0.25,-0.75,0,0,0,0,0,-0.75,0,0,0,0,0,-0.25,0,0.25,0.25,0.5,0.75,0,0,0,1],
            [-0.25,-0.5,0,0,0,0,0,0,0,0,1,0,-0.75,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.75,1]
          ]
        )
        
    # 2×1（投資家のビュー）
    Q = np.array([[0.05], [0.07]])
    
    # 2×1（投資家のビューへの確信度）
    omega = np.array(
              [
                [0.001065383, 0],
                [0, 0.001851738]
              ]
            )
            
    for i in daily.columns:
      daily[i] = cp[i].pct_change()
    
    sigma = daily[:245].cov().as_matrix()
    
    delta = 25
    tau = 0.05
    ganma = 0.3
    
    count = 0
    w_matrix =[]
    
    for index, row in weight[245:].iterrows():
      row_list = list(row)
      count += 1
      w_out = []
      for num, i in enumerate(row_list):
        w_in = []
        w_in.append(row_list[num])
        w_out.append(w_in)
      w = np.array(w_out)
      
      r_eq = delta * np.dot(sigma, w)
      r_posterior = r_eq + np.dot(np.dot(tau*np.dot(sigma, P.T), np.linalg.inv(tau*np.dot(np.dot(P, sigma), P.T) + omega)), (Q - np.dot(P, r_eq)))
      sigma_posterior = sigma + tau * sigma - tau * np.dot(np.dot(np.dot(sigma, P.T), np.linalg.inv(tau*np.dot(np.dot(P, sigma), P.T)+omega)), tau*np.dot(P, sigma))
      w_posterior = np.dot(np.linalg.inv(delta * sigma_posterior), r_posterior).T
      w_matrix = np.append(w_matrix, w_posterior)
      
    w_matrix_rt = w_matrix.reshape(count, 30).T
    
    for sym, val in zip(weight[245:].keys(), w_matrix_rt):
      post = pd.Series(val, index=cp[245:], name = sym)
      df_w_posterior[sym] = post.values
      
    buy_sig = cp[(cp > 0)]
    sell_sig = cp[(cp < 0)]
    
    buy_sig[~pd.isnull(buy_sig)] = 1.0
    sell_sig[~pd.isnull(buy_sig)] = 0.0
    sell_sig[~pd.isnull(sell_sig)] = 1.0
    
    return {
      "buy:sig": buy_sig,
      "sell:sig": sell_sig,
      "weight:g2": weight,
      "weight_post:g2": df_w_posterior,
      "weight_comp:g2": (1 - ganma) * weight + ganma * df_w_posterior,
      "weight_diff:g2": weight - df_w_posterior,
    }
    
  ctx.regist_signal("my_signal", _my_signal)
  
def handle_signals(ctx, date, current):
  df = current.copy()
  buy = df["buy:sig"].dropna()
  buy = buy[buy != 0]
  
  if not buy.empty:
    for (sym, val) in buy.items():
      sec = ctx.getSecurity(sym)
      sec.order_target_percent(df["weight_comp:g2"][sym])
```
  </div>
</details>

```python
def _my_signal(data):
  #終値(株式分割調整後)
  cp = data["close_price_adj"].fillna(method="ffill")
  #時価総額
  mkt_cap = data["market_capitalization"].fillna(method="ffill")
  #時価総額の合計
  sum_mkt_cap = pd.Series(data=0, index=cp.index)
  #時価総額ベースの重み
  weight = pd.DataFrame(data=0, index=cp.index, columns=[])
  #日次リターン（終値の差）
  daily = pd.DataFrame(data=0, index=cp.index, columns=cp.columns)
  #調整後の重み
  df_w_posterior = pd.DataFrame(data=0, index=cp.index[245:], columns=cp.columns)
```

色々な変数の容れ物を用意しています．

```python
  for sym, val in cp.items():
    sum_mkt_cap += mkt_cap[sym]
```

各銘柄の時価総額を足していくことで，全体の時価総額の合計を出しています．

```python
  for sym, val in cp.items():
    weight[sym] = mkt_cap[sym] / sum_mkt_cap
```

各銘柄の時価総額を全体の時価総額の合計で割ることで時価総額ベースの重みを出しています．

```python
# 2x30 (投資家のビュー，相対的な影響度)
P = np.array(
      [
        [-0.5,-0.25,0,0,0,0,0,-0.25,-0.75,0,0,0,0,0,-0.75,0,0,0,0,0,-0.25,0,0.25,0.25,0.5,0.75,0,0,0,1],
        [-0.25,-0.5,0,0,0,0,0,0,0,0,1,0,-0.75,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.75,1]
      ]
    )

# 2×1（投資家のビュー，影響度の大きさ）
Q = np.array([[0.05], [0.07]])

# 2×1（投資家のビューへの確信度）
omega = np.array(
          [
            [0.001065383, 0],
            [0, 0.001851738]
          ]
        )
```

投資家のビューに適当な見通しを与えています．

```python
for i in daily.columns:
  #変化の比率を計算する
  daily[i] = cp[i].pct_change()
  #リターンの共分散行列  
  sigma = daily[:245].cov().as_matrix()

  #投資家のリスク回避度    
  delta = 25
  #共分散行列に対する信頼区間
  tau = 0.05
  #調整後の重みの適用割合
  ganma = 0.3
  
  #テスト期間の営業日の数をカウントする
  count = 0
  #計算途中の重みの容れ物
  w_matrix =[]
```

何種類かのパラメーターに適当な値を与えています．

```python
for index, row in weight[245:].iterrows():
  row_list = list(row)
  count += 1
  #外側のリスト
  w_out = []

  for num, i in enumerate(row_list):
    #内側のリスト
    w_in = []
    w_in.append(row_list[num])
    w_out.append(w_in)
  #調整前重み
  w = np.array(w_out)
```
続く計算に渡すのに適した形に変形しています．

以降の数式で登場する記号は以下のとおりです．（N は銘柄の数，K は投資家のビューの数）

$\Pi$ 均衡期待リターン（N 次元ベクトル）
$\Pi'$ 投資家のビューを考慮して更新された期待リターン（N 次元ベクトル）
$w$ 当初ウェイト（N 次元ベクトル）
$w'$ 事後ウェイト（N 次元ベクトル）
$\Sigma$ リターンの共分散行列（N x N 行列）
$P$ 投資家のビューの相対的な影響度を示す行列（K x N 行列）
$Q$ 投資家のビューの影響度の大きさを示す行列（K 次元ベクトル）
$\Omega$ 投資家のビューに対する確信度を示す行列（K x K 行列）
$\delta$ パラメータ（投資家のリスク回避度を表す）
$\tau$ パラメータ（共分散行列に対する信頼区間を表す）

```python
  #均衡期待リターン
  r_eq = delta * np.dot(sigma, w)
```

```math
\Pi = \delta\Sigma w
```

```python
  #均衡期待リターンに投資家のビューを反映する
  r_posterior = r_eq + np.dot(np.dot(tau*np.dot(sigma, P.T), np.linalg.inv(tau*np.dot(np.dot(P, sigma), P.T) + omega)), (Q-np.dot(P, r_eq)))
```

```math
\Pi' = \Pi + \tau \Sigma P^T \left( P \tau \Sigma P^T + \Omega \right)^{-1} \left( Q - P\Pi\right)
```

```python
  #リターンの共分散行列に投資家のビューを反映する
  sigma_posterior = sigma + tau * sigma - tau * np.dot(np.dot(np.dot(sigma, P.T), np.linalg.inv(tau*np.dot(np.dot(P, sigma), P.T)+omega)), tau*np.dot(P, sigma))
```

```math
\Sigma' = \Sigma + \tau\Sigma - \tau\Sigma P^T \left( P \tau \Sigma P^T + \Omega \right)^{-1} P\tau\Sigma
```

```python
  #Forward Optimizationをして最適重みを求める
  w_posterior = np.dot(np.linalg.inv(delta * sigma_posterior), r_posterior).T
```

```math
w' = \Pi' (\delta\Sigma')^{-1}
```

```python
  #計算途中の重みの容れ物に加えていく
  w_matrix = np.append(w_matrix, w_posterior)
```

```python
#調整後重みを変形する
w_matrix_rt = w_matrix.reshape(count, 30).T
    
for sym, val in zip(weight[245:].keys(), w_matrix_rt):
  post = pd.Series(val, index=cp[245:], name = sym)
  #pandas.DataFrameに収める
  df_w_posterior[sym] = post.values
```

調整前重みのデータの形式と調整後重みのデータの形式を同じものとするために，このような操作をしています．

```python
#適当な売買シグナル
buy_sig = cp[(cp > 0)]
sell_sig = cp[(cp < 0)]
    
buy_sig[~pd.isnull(buy_sig)] = 1.0
sell_sig[~pd.isnull(buy_sig)] = 0.0
sell_sig[~pd.isnull(sell_sig)] = 1.0
```

```python
return {
  #買いシグナル
  "buy:sig": buy_sig,
  #売りシグナル
  "sell:sig": sell_sig,
  #調整前重み
  "weight:g2": weight,
  #調整後重み
  "weight_post:g2": df_w_posterior,
  #調整前重みと調整後重みの混合重み
  "weight_comp:g2": (1 - ganma) * weight + ganma * df_w_posterior,
  #調整前重みと調整後重みの差
  "weight_diff:g2": weight - df_w_posterior,
}
```

テスト実行後に個別銘柄詳細の欄で，重みの変化などを確認することができます．

```python
def handle_signals(ctx, date, current):
  df = current.copy()
  #買いシグナルの容れ物
  buy = df["buy:sig"].dropna()
  buy = buy[buy != 0]
  
  if not buy.empty:
    for (sym, val) in buy.items():
      #銘柄の特定
      sec = ctx.getSecurity(sym)
      #特定銘柄の総保有額が総資産評価額に対して指定の割合となるように注文をする
      sec.order_target_percent(df["weight_comp:g2"][sym])
```

ここでは，買いシグナルに沿って注文を実行しています．

# おわりに
下記のURLからQuantX Factory上で今回紹介したコードのバックテストを行うことが出来ます．
https://factory.quantx.io/developer/9deac3f08ff44b4ca49be39f1888578f

# 参考文献
[1] Fischer Black and Robert Litterman, "Global Portfolio Optimization", Financial Analysts Journal, Vol. 48, No. 5, pp. 28-43, 1992.

[2] Guangliang He and Robert Litterman, "The Intuition Behind Black-Litterman Model Portfolios", Goldman Sachs Investment Management Research, 1999.

[3] Stephen Satchell and Alan Scowcroft, "A demystification of the Black–Litterman model: Managing quantitative and traditional portfolio construction", Journal of Asset Management, Vol. 1, No. 1, pp. 138-150, 2000.

[4] Thomas Idzorek, "A step-by-step guide to the Black-Litterman model", Forecasting Expected Returns in the Financial Markets, pp. 17–38, 2007.

[5] Harry Markowitz, "Portfolio Selection", The Journal of Finance, Vol. 7, No. 1, pp. 77-91, 1952.

[6] William F. Sharpe, "Capital Asset Prices: A Theory of Market Equilibrium under Conditions of Risk", The Journal of Finance, Vol. 19, No. 3, pp. 425-442, 1964.

また，以下の記事を参考にしています．
[ブラック・リッターマンモデルによる資産配分を解説してみる（Pythonによる実行例つき）](https://qiita.com/nokomitch/items/0d1812763114e6266bf3)

#免責注意事項
この記事のコードや知識を使った実際の取引で生じた損益に関しては一切の責任を負いかねますので御了承下さい．

[^1]: https://www.nli-research.co.jp/files/topics/37745_ext_18_0.pdf
[^2]: https://en.wikipedia.org/wiki/Black%E2%80%93Litterman_model