【証券分析10-3】リスクパリティポートフォリオを作ってみる

前回は、リスクパリティポートフォリオを作るにあたって、どのような条件が満たされていればよいのか、どんな問題を解けば解が得られるのか見ていきました。そこで今回は、前回、明らかにした最適化問題をPythonで実際に解き、シミュレーションを行っていきたいと思います。

解くべき最適化問題

リスクパリティポートフォリオをためには、以下の最適化問題を解く必要があることを前回は見ていきました。
ここからは、Pythonでこの最適化問題を解いていきます。

$$\text{minimize } \sum_{i=1}^{n}(RC_{i}-\frac{\sigma_{p}}{n})^{2}$$

$$\text{subject to } w_{1}≧0$$

$$資産iのリスク寄与=RC_{i} = w_{i}\frac{\partial\sigma_{p}}{\partial w_{i}}$$

$$=\frac{w_{i}\sum_{j=1}^{n}w_{j}\sigma_{ij}}{\sigma_{p}}$$

リスクパリティポートフォリオ導出

まず、素直に上記の最適化問題を関数にした場合、以下のような関数を作成することができます。
今回は、Pythonのscipyというライブラリを使って最適化を行います。
コードの詳細な解説は、記事の下に記載しています。
実際のコードは以下のGoogleColaboratoryでも確認できます。
コードはこちら

import scipy.optimize as op
import pandas as pd
import numpy as np      

def cal_risk_parity(Sigma):
  #資産数
  n = Sigma.shape[0]
  #ポートフォリオの分散を計算する関数
  def calculate_portfolio_var(w,Sigma):
    w = np.matrix(w)
    return np.dot(np.dot(w,Sigma),w.T)

  #RCの計算
  def calculate_risk_contribution(w,Sigma):
    w = np.matrix(w)
    sigma_p = np.sqrt(calculate_portfolio_var(w,Sigma))
    RC = np.multiply(np.dot(Sigma,w.T),w.T)/sigma_p
    return RC

  #目的関数の設定
  def risk_parity_objective(w):
    sigma_p =  np.sqrt(calculate_portfolio_var(w,Sigma))
    risk_target = np.asmatrix(np.multiply(sigma_p,[1./n]*n))
    RC = calculate_risk_contribution(w,Sigma)
    obj = sum(np.square(RC-risk_target.T))
    return obj

  #制約条件
  cons = [{'type': 'eq', 'fun': lambda w: np.sum(w) - 1}] #ウェイト合計は1
  bnds = [(0, None)] * n  #各資産のウェイトは0以上
  #初期値
  x0 = [1./n] * n
  #最適化実行
  opts = op.minimize(risk_parity_objective, x0=x0, method='SLSQP', bounds=bnds, constraints=cons)
  return opts['x']

2資産の場合

これを2資産の場合で実行してみると、

#資産1のリスク13%,3のリスク6%、相関係数0.5の場合
s1 = 0.13
rho = 0.5
s2 = 0.06
Sigma = np.array([[s1**2, rho*s1*s2],[rho*s1*s2, s2**2]])
cal_risk_parity(Sigma)

array([0.31578841, 0.68421159])

となり、これは2資産のリスクから求める結果と同じになります。

#2資産の場合は簡単に解が得られる
print(s2/(s1+s2))
print(s1/(s1+s2))

0.3157894736842105
0.6842105263157895

RCの値などを出力してみると、

ポートフォリオのリスク:[[0.07111077]]
ターゲットリスク:[[0.03555538 0.03555538]]
RC: [[0.03556537]
 [0.03554539]]
目的関数；[[1.99577376e-10]]

このように最終的には、両方の資産のRCがほぼ同じになっていることを確認できます。

ちなみに、50%ずつの場合は、以下のようになります。

ポートフォリオのリスク:[[0.08411302]]
ターゲットリスク:[[0.04205651 0.04205651]]
RC: [[0.06182158]
 [0.02229144]]
目的関数；[[0.00078132]]

4資産の場合

次に　4資産の場合を見てみましょう。

#4資産の場合
s1 = 0.13
s2 = 0.06
s3 = 0.03
s4 = 0.08
rho12 = 0.5
rho13 = 0.2
rho14 = -0.3
rho23 = 0.4
rho24 = -0.1
rho34 = 0.2
Sigma = np.array([
    [s1**2, rho12*s1*s2, rho13*s1*s3, rho14*s1*s4],
     [rho12*s1*s2, s2**2, rho23*s2*s3, rho24*s2*s4],
    [rho13*s1*s3, rho23*s2*s3, s3*s3, rho34*s3*s4],
    [rho14*s1*s4, rho24*s2*s4, rho34*s3*s4, s4*s4]
    ])
w = cal_risk_parity(Sigma)
w

array([0.15236437, 0.26592141, 0.27711471, 0.30459951])

この時のRCは、

n = Sigma.shape[0]
sigma_p =  np.sqrt(calculate_portfolio_var(w,Sigma))
print(f'ポートフォリオのリスク:{sigma_p}')
risk_target = np.asmatrix(np.multiply(sigma_p,[1./n]*n))
print(f'ターゲットリスク:{risk_target}')
RC = calculate_risk_contribution(w,Sigma)
print(f'RC: {RC}')
obj = sum(np.square(RC-risk_target.T))
print(f'目的関数；{obj}')

ポートフォリオのリスク:[[0.03887803]]
ターゲットリスク:[[0.00971951 0.00971951 0.00971951 0.00971951]]
RC: [[0.0112784 ]
 [0.01097701]
 [0.00503164]
 [0.01159098]]
目的関数；[[0.2948993]]

このようになっていて、あまり近くありません。
思ったよりもずれが大きくなってしまうが、これはすでに目的関数が十分に小さくなっているとして、計算が終了していることに問題があります。
そこで、目的関数を％単位に変換sh知恵から誤差を算出するように変更して実行してみます。

改良版リスクパリティ関数

import scipy.optimize as op
import pandas as pd
import numpy as np      

def cal_risk_parity2(Sigma):
  #資産数
  n = Sigma.shape[0]
  #ポートフォリオの分散を計算する関数
  def calculate_portfolio_var(w,Sigma):
    w = np.matrix(w)
    return np.dot(np.dot(w,Sigma),w.T)

  #RCの計算
  def calculate_risk_contribution(w,Sigma):
    w = np.matrix(w)
    sigma_p = np.sqrt(calculate_portfolio_var(w,Sigma))
    RC = np.multiply(np.dot(Sigma,w.T),w.T)/sigma_p
    return RC

  #目的関数の設定
  def risk_parity_objective(w):
    sigma_p =  np.sqrt(calculate_portfolio_var(w,Sigma))
    risk_target = np.asmatrix(np.multiply(sigma_p,[1./n]*n))
    RC = calculate_risk_contribution(w,Sigma)
    obj = sum(np.square(RC*100-risk_target.T*100)) #sum(np.square(RC-risk_target.T))も同じことだが、資産が増えるとそれだけで誤差が小さくなるので、100倍しておく
    return obj

  #制約条件
  cons = [{'type': 'eq', 'fun': lambda w: np.sum(w) - 1}] #ウェイト合計は1
  bnds = [(0, None)] * n  #各資産のウェイトは0以上
  #初期値
  x0 = [1./n] * n
  #最適化実行
  opts = op.minimize(risk_parity_objective, x0=x0, method='SLSQP', bounds=bnds, constraints=cons)
  return opts['x']

w = cal_risk_parity2(Sigma)
print(w)
n = Sigma.shape[0]
sigma_p =  np.sqrt(calculate_portfolio_var(w,Sigma))
print(f'ポートフォリオのリスク:{sigma_p}')
risk_target = np.asmatrix(np.multiply(sigma_p,[1./n]*n))
print(f'ターゲットリスク:{risk_target}')
RC = calculate_risk_contribution(w,Sigma)
print(f'RC: {RC}')
obj = sum(np.square(RC*100-risk_target.T*100))
print(f'目的関数；{obj}')

[0.12406626 0.21518608 0.41671485 0.2440328 ]
ポートフォリオのリスク:[[0.03522153]]
ターゲットリスク:[[0.00880538 0.00880538 0.00880538 0.00880538]]
RC: [[0.00880473]
 [0.00880641]
 [0.00880109]
 [0.0088093 ]]
目的関数；[[3.52892111e-07]]

このようにかなりリスク寄与が均一になりました。

分散ベースで考える場合

前回は、偏微分を用いた数式だけでなく、分散をベースにした考え方も示しました。
そちらの場合は、以下のような実装になります。
基本的には同じですが、RCの計算をするときに分散の状態のままで計算する点と目的関数の設定で、分散を使っている部分が違いになります。

import scipy.optimize as op
import pandas as pd
import numpy as np      

def cal_risk_parity3(Sigma):
  #資産数
  n = Sigma.shape[0]
  #ポートフォリオの分散を計算する関数
  def calculate_portfolio_var(w,Sigma):
    w = np.matrix(w)
    return np.dot(np.dot(w,Sigma),w.T)

  #RCの計算
  def calculate_risk_contribution(w,Sigma):
    w = np.matrix(w)
    var_p = calculate_portfolio_var(w,Sigma)
    RC = np.multiply(np.dot(Sigma,w.T),w.T)
    return RC

  #目的関数の設定
  def risk_parity_objective(w):
    var_p =  calculate_portfolio_var(w,Sigma)
    risk_target = np.asmatrix(np.multiply(var_p,[1./n]*n))
    RC = calculate_risk_contribution(w,Sigma)
    obj = sum(np.square(RC*100-risk_target.T*100)) #sum(np.square(RC-risk_target.T))も同じことだが、資産が増えるとそれだけで誤差が小さくなるので、100倍しておく
    return obj

  #制約条件
  cons = [{'type': 'eq', 'fun': lambda w: np.sum(w) - 1}] #ウェイト合計は1
  bnds = [(0, None)] * n  #各資産のウェイトは0以上
  #初期値
  x0 = [1./n] * n
  #最適化実行
  opts = op.minimize(risk_parity_objective, x0=x0, method='SLSQP', bounds=bnds, constraints=cons)
  return opts['x']

この場合でも同様の結果を得られます。

w = cal_risk_parity2(Sigma)
print(w)

[0.12406626 0.21518608 0.41671485 0.2440328 ]

このようにscipyを使うことで、簡単にリスクパリティウェイトを求めることができます。

コードの解説

必要なライブラリのインポート

import scipy.optimize as op
import pandas as pd
import numpy as np

まず、最適化のために必要なライブラリを導入します。

ポートフォリオ最適化関数の定義

def cal_risk_parity2(Sigma):

ポートフォリオ最適化関数を定義します。引数として共分散行列（Sigma）を取ります。

資産数の取得

n = Sigma.shape[0]

共分散行列の行数から資産数を取得します。

ポートフォリオ分散計算関数の定義

def calculate_portfolio_var(w,Sigma):
    w = np.matrix(w)
    return np.dot(np.dot(w,Sigma),w.T)

関数 calculate_risk_contribution は、リスク寄与を計算するものです。引数としてはウェイトベクトル（w）と共分散行列（Sigma）です。
ウェイトベクトル w を行列に変換して、標準偏差を求めます。

リスク寄与計算関数の定義

def calculate_risk_contribution(w, Sigma):

リスク寄与を計算する関数を定義します。ウェイトベクトル（w）と共分散行列（Sigma）を引数に取ります。

sigma_p = np.sqrt(calculate_portfolio_var(w, Sigma))

calculate_portfolio_var 関数を使用してポートフォリオの分散を計算し、その平方根を取得します。

RC = np.multiply(np.dot(Sigma, w.T), w.T) / sigma_p

共分散行列 Sigma とウェイトベクトル w を用いてリスク寄与を計算します。計算結果をポートフォリオの標準偏差で割ります。

目的関数の設定

def risk_parity_objective(w):

目的関数 risk_parity_objective は、リスク平等ポートフォリオの構築において最適なウェイトを見つけるための関数です。引数としてはウェイトベクトル（w）が渡されます。

sigma_p = np.sqrt(calculate_portfolio_var(w, Sigma))

calculate_portfolio_var 関数を使用してポートフォリオの分散を計算し、その平方根を取得します。

risk_target = np.asmatrix(np.multiply(sigma_p, [1. / n] * n))

目標リスクを計算し、行列に変換します。目標リスクは資産数に基づいてポートフォリオの平均リスクを均等に分配するためのものです。

RC = calculate_risk_contribution(w, Sigma)

先ほど定義した calculate_risk_contribution 関数を使用してリスク寄与を計算します。

obj = sum(np.square(RC * 100 - risk_target.T * 100))

リスク寄与と目標リスクの誤差を算出し、その二乗の合計を目的関数とします。誤差を100倍しているのは、誤差が小さくなりすぎることを防ぐためです。

制約条件の設定

cons = [{'type': 'eq', 'fun': lambda w: np.sum(w) - 1}]
bnds = [(0, None)] * n

ウェイトの合計が1である制約条件と各資産のウェイトが0以上である制約条件を設定します。

最適化の実行

opts = op.minimize(risk_parity_objective, x0=x0, method='SLSQP', bounds=bnds, constraints=cons)
return opts['x']

最適化を実行し、最適なウェイトを取得します。