このページは、nAGライブラリのJupyterノートブックExampleの日本語翻訳版です。オリジナルのノートブックはインタラクティブに操作することができます。
nAGライブラリを使用したポートフォリオ最適化における堅牢な線形計画法
はじめに
損失リスク制約付きの古典的なポートフォリオ問題を考えます： \[\begin{equation}\label{prob} \begin{array}{ll} \underset{x\in\Re^n}{\mbox{maximize}} & \bar{r}^Tx\\[0.6ex] \mbox{subject to} & \mbox{Prob}(r^Tx\leq\alpha)\leq\beta,\\[0.6ex] & \sum_{i=1}^n x_i=1,\\[0.6ex] & x\geq0. \end{array} \end{equation}\] ここで、資産の収益率 \(r\in\Re^n\) は平均 \(\bar{r}\in\Re^n\) と共分散 \(V\in\Re^{n\times n}\) を持つガウス分布に従うと仮定し、\(\alpha\) は与えられた望ましくない収益レベル（例：過度の損失）、\(\beta\) は与えられた最大確率です。
問題 (prob) は以下の形式の二次錐計画問題（SOCP）に変換でき、nAGライブラリのMark 27で導入されたSOCPソルバーを使用して簡単に解けることを示します： \[\begin{equation}\label{SOCP} \begin{array}{ll} \underset{x\in\Re^n}{\mbox{minimize}} & c^Tx\\[0.6ex] \mbox{subject to} & l_A\leq Ax\leq u_A,\\[0.6ex] & l_x\leq x\leq u_x,\\[0.6ex] & x\in{\cal K}, \end{array} \end{equation}\] ここで \(A\in\Re^{m\times n}\), \(l_A, u_A\in\Re^m\), \(c, l_x, u_x\in\Re^n\) は問題データであり、\(\cal K={\cal K}^{n_1}\times\cdots\times{\cal K}^{n_r}\times\Re^{n_l}\) で \({\cal K}^{n_i}\) は以下のように定義される二次錐または回転二次錐のいずれかです：確率制約のより詳細な検討
\(u=r^Tx\)とし、\(\sigma=x^TVx\)をその分散とすると、問題(prob)の確率制約は以下のように書くことができます \[ \mbox{Prob}\left(\frac{u-\bar{u}}{\sqrt{\sigma}}\leq\frac{\alpha-\bar{u}}{\sqrt{\sigma}}\right)\leq\beta. \] \((u-\bar{u})/\sqrt{\sigma}\)は標準正規分布に従う確率変数であることに注意してください。上記の確率は単に\(\Phi((\alpha-\bar{u})/\sqrt{\sigma})\)です。ここで、 \[ \Phi(z) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^z e^{-t^2/2}dt \] は標準正規分布の累積分布関数です。したがって、問題(prob)の確率制約は以下のように書くことができます \[ \frac{\alpha-\bar{u}}{\sqrt{\sigma}}\leq\Phi^{-1}(\beta) \] または、同等に、 \[\begin{equation}\label{label_1} \bar{u} + \Phi^{-1}(\beta)\sqrt{\sigma}\geq\alpha. \end{equation}\] \(\bar{u}=\bar{r}^Tx\)と\(\sigma=x^TVx=\|Fx\|^2\)（\(V=F^TF\)と因数分解して）から、(label_1)は以下と同等です \[\begin{equation}\label{label_2} \bar{r}^Tx + \Phi^{-1}(\beta)\|Fx\|\geq\alpha. \end{equation}\] \(\beta\)の値に応じて、(label_2)は凸または凹になる可能性があります。\(\beta\leq0.5\)と設定すると（これはリスク管理において合理的です）、(label_2)は凸となり、以下のように書くことができます \[\begin{equation}\label{label_3} \bar{r}^Tx + \Phi^{-1}(\beta)t=\alpha ~ \mbox{and} ~ \|Fx\|\leq t. \end{equation}\] \(Fx=y\)とすると、(label_3)の\(\|Fx\|=\|y\|\leq t\)は二次錐制約(SOC)に完全に適合することに注意してください。したがって、二次錐計画法で解くことができます。問題(prob)の最終的な等価なSOCPモデルは以下のようになります \[\begin{equation}\label{final} \begin{array}{ll} \underset{x\in\Re^n}{\mbox{maximize}} & \bar{r}^Tx\\[0.6ex] \mbox{subject to} & \sum_{i=1}^n x_i=1,\\[0.6ex] & Fx = y,\\[0.6ex] & \bar{r}^Tx + \Phi^{-1}(\beta)t=\alpha, \\[0.6ex] & (t;y)\in{\cal K}^{n+1}_q,\\[0.6ex] & x\geq0. \end{array} \end{equation}\]
nAGライブラリの使用
nAG SOCPソルバーを使用して問題(prob)をモデル化し、等価なSOCP問題(final)を解く方法を示します。
# ユーティリティライブラリとnAGライブラリをインポートする
import numpy as np
import math as mt
from scipy.stats import norm
from naginterfaces.library import opt
from naginterfaces.library import lapackeig例として、\(\alpha=0.0001\)、\(\beta=0.05\)とし、\(8\)個の資産に対して\(\bar{r}\)と\(V\)をランダムに生成します。以下のようになります。
# アルファとベータを設定する
alpha = 0.0001
beta = 0.05
# 乱数シードを固定する
np.random.seed(9)
# 資産数
n_assets = 8
# 期待リターンのベクトル
r = np.ones(n_assets)*.02 + np.random.rand(n_assets)*.15
# 共分散行列
V = np.matrix(np.random.randn(n_assets, n_assets))
V = V.T * V
V = V / np.max(np.abs(np.diag(V))) * .2nAGライブラリの固有値分解を使用して、\(V=F^TF\)を因数分解します。
# 注意: Vがスパース行列として入力される場合、スパース因子分解を使用することができます
U, lamda = lapackeig.dsyevd('V', 'L', V)
# 正の固有値と対応する固有ベクトルを求める
i = 0
k = 0
F = []
while i<len(lamda):
if lamda[i] > 0:
F = np.append(F, mt.sqrt(lamda[i])*U[:,i])
k += 1
i += 1
F = F.reshape((k, n_assets))モデリングのために、nAG SOCP ソルバーは目的関数と制約条件のためのいくつかの入力引数を必要とします。ここで、nAG SOCP ソルバーに供給するデータを初期化します。
# 変数の数
n = n_assets
# 制約条件の数
m = 0
# 目的関数の係数 c
c = np.full(n, 0.0, float)
# 変数の境界
blx = np.full(n, -1.e20, float)
bux = np.full(n, 1.e20, float)
# 線形制約 bu <= Ax <= bu
# 座標リスト形式(COO)のA
irowa = np.empty(0, int)
icola = np.empty(0, int)
a = np.empty(0, float)
# Axの境界
bl = np.empty(0, float)
bu = np.empty(0, float)
# 円錐制約
ctype = []
group = []プロセス中に補助変数と制約を追加するため、モデル内の変数と制約の数を追跡し、最新の問題サイズを維持する必要があります。
# 最新の問題サイズを初期化する
n_up = n
m_up = mここで、目的関数と制約条件を1つずつ追加する際に、上記のデータを継続的に修正していきます。最初は目的関数の係数 \(c\) です。
# 目的関数を追加 min -r'x
c = -rそして、ロングオンリー制約を追加します。
# 線形制約の数が1増加します
m += 1
# x の下限を0に設定
blx[0:n] = np.zeros(n)
# 和を追加 sum(x) = 1
irowa = np.append(irowa, np.full(n, m_up+1, dtype=int))
icola = np.append(icola, np.arange(1, n+1))
a = np.append(a, np.full(n, 1.0, dtype=float))
bl = np.append(bl, np.full(1, 1.0, dtype=float))
bu = np.append(bu, np.full(1, 1.0, dtype=float))以下に確率制約を追加します \[ Fx = y,~ \bar{r}^Tx + \Phi^{-1}(\beta)t=\alpha ~及び~(t;y)\in{\cal K}^{n+1}_q. \]
# ベータ分布の分位関数を取得
quantile = norm.ppf(beta)
# 最新の問題サイズ
m_up = m
n_up = n
# その後、k + 1個の変数を以下と共に加える必要があります
# k + 1個の線形制約と1つの二次錐制約
# モデルを拡大する
n = n + k + 1
m = m + k + 1
# 追加されたすべての補助変数は目的関数に関与しない
c = np.append(c, np.zeros(k+1))
# xの範囲を拡大し、新たに追加されたk+1個の変数に無限大の範囲を追加する
blx = np.append(blx, np.full(k+1, -1.e20, dtype=float))
bux = np.append(bux, np.full(k+1, 1.e20, dtype=float))
# 線形制約の拡大
# F (COO)の疎パターン
row, col = np.nonzero(F)
val = F[row, col]
# 1ベースに変換し、行を m だけ下に移動する
# Fx = y と Phi^-1(beta)*t + r_bar'x - alpha = 0 を追加
# [x,t,y]
row = row + 1 + m_up
col = col + 1
row = np.append(row, np.arange(m_up+1, m_up+k+1+1))
col = np.append(np.append(col, np.arange(n_up+2, n_up+k+1+1)), n_up+1)
val = np.append(val, np.append(np.full(k, -1.0, dtype=float), quantile))
irowa = np.append(irowa, row)
icola = np.append(icola, col)
a = np.append(a, val)
bl = np.append(bl, np.append(np.zeros(k), alpha))
bu = np.append(bu, np.append(np.zeros(k), alpha))
# xの係数は Phi^-1(beta)*t + r'x - alpha = 0 において
irowa = np.append(irowa, np.full(n_assets, m_up+k+1, dtype=int))
icola = np.append(icola, np.arange(1, n_assets+1))
a = np.append(a, r)
# 円錐制約の拡大
ctype.extend('Q')
group_temp = np.arange(n_up+1, n_up+k+1+1)
group.append(group_temp)これで、必要なすべてのデータが準備できました。それらをnAG SOCP ソルバーに入力して解を求めます。
# 問題ハンドルを作成する
handle = opt.handle_init(n)
# 目的関数を設定する
opt.handle_set_linobj(handle, c)
# 箱制約を設定する
opt.handle_set_simplebounds(handle, blx, bux)
# 線形制約を設定する
opt.handle_set_linconstr(handle, bl, bu, irowa, icola, a)
# 円錐制約を設定する
i = 0
while i<len(ctype):
opt.handle_set_group(handle, ctype[i], 0, group[i])
i += 1
# オプションを設定
for option in [
'Print Options = NO',
'Print File = 1',
'SOCP Scaling = A'
]:
opt.handle_opt_set(handle, option)
# 内点法SOCPソルバーを呼び出す
slt = opt.handle_solve_socp_ipm(handle)naginterfaces.base.opt.handle_solve_socp_ipm:
naginterfaces.base.opt.handle_solve_socp_ipm: ------------------------------------------------
naginterfaces.base.opt.handle_solve_socp_ipm: E04PT, Interior point method for SOCP problems
naginterfaces.base.opt.handle_solve_socp_ipm: ------------------------------------------------
naginterfaces.base.opt.handle_solve_socp_ipm:
naginterfaces.base.opt.handle_solve_socp_ipm: Problem Statistics
naginterfaces.base.opt.handle_solve_socp_ipm: No of variables 17
naginterfaces.base.opt.handle_solve_socp_ipm: free (unconstrained) 9
naginterfaces.base.opt.handle_solve_socp_ipm: bounded 8
naginterfaces.base.opt.handle_solve_socp_ipm: No of lin. constraints 10
naginterfaces.base.opt.handle_solve_socp_ipm: nonzeroes 89
naginterfaces.base.opt.handle_solve_socp_ipm: No of quad.constraints 0
naginterfaces.base.opt.handle_solve_socp_ipm: No of cones 1
naginterfaces.base.opt.handle_solve_socp_ipm: biggest cone size 9
naginterfaces.base.opt.handle_solve_socp_ipm: Objective function Linear
naginterfaces.base.opt.handle_solve_socp_ipm:
naginterfaces.base.opt.handle_solve_socp_ipm: Presolved Problem Measures
naginterfaces.base.opt.handle_solve_socp_ipm: No of variables 17
naginterfaces.base.opt.handle_solve_socp_ipm: No of lin. constraints 10
naginterfaces.base.opt.handle_solve_socp_ipm: nonzeroes 89
naginterfaces.base.opt.handle_solve_socp_ipm: No of cones 1
naginterfaces.base.opt.handle_solve_socp_ipm:
naginterfaces.base.opt.handle_solve_socp_ipm:
naginterfaces.base.opt.handle_solve_socp_ipm: ------------------------------------------------------------------------
naginterfaces.base.opt.handle_solve_socp_ipm: it| pobj | dobj | p.inf | d.inf | d.gap | tau | I
naginterfaces.base.opt.handle_solve_socp_ipm: ------------------------------------------------------------------------
naginterfaces.base.opt.handle_solve_socp_ipm: 0 -1.13221E+00 0.00000E+00 8.43E-01 2.14E-01 4.45E+00 1.0E+00
naginterfaces.base.opt.handle_solve_socp_ipm: 1 -2.82490E-01 -8.98481E-02 1.91E-01 4.87E-02 1.01E+00 1.3E+00
naginterfaces.base.opt.handle_solve_socp_ipm: 2 -1.32896E-01 -9.68963E-02 3.73E-02 9.50E-03 1.97E-01 1.4E+00
naginterfaces.base.opt.handle_solve_socp_ipm: 3 -1.10308E-01 -9.31071E-02 1.83E-02 4.66E-03 9.65E-02 1.4E+00
naginterfaces.base.opt.handle_solve_socp_ipm: 4 -1.02701E-01 -8.92455E-02 1.10E-02 2.81E-03 5.83E-02 1.1E+00
naginterfaces.base.opt.handle_solve_socp_ipm: 5 -8.78903E-02 -8.59020E-02 1.62E-03 4.12E-04 8.55E-03 1.1E+00
naginterfaces.base.opt.handle_solve_socp_ipm: 6 -8.58224E-02 -8.53194E-02 3.89E-04 9.90E-05 2.05E-03 1.0E+00
naginterfaces.base.opt.handle_solve_socp_ipm: 7 -8.54027E-02 -8.51613E-02 1.71E-04 4.34E-05 9.01E-04 9.6E-01
naginterfaces.base.opt.handle_solve_socp_ipm: 8 -8.50652E-02 -8.50604E-02 3.51E-06 8.93E-07 1.85E-05 9.7E-01
naginterfaces.base.opt.handle_solve_socp_ipm: 9 -8.50588E-02 -8.50583E-02 3.86E-07 9.84E-08 2.04E-06 9.7E-01
naginterfaces.base.opt.handle_solve_socp_ipm: 10 -8.50580E-02 -8.50580E-02 9.88E-09 2.51E-09 5.21E-08 9.7E-01
naginterfaces.base.opt.handle_solve_socp_ipm: 11 -8.50580E-02 -8.50580E-02 3.19E-10 1.47E-10 8.04E-10 9.7E-01
naginterfaces.base.opt.handle_solve_socp_ipm: ------------------------------------------------------------------------------
naginterfaces.base.opt.handle_solve_socp_ipm: Status: converged, an optimal solution found
naginterfaces.base.opt.handle_solve_socp_ipm: ------------------------------------------------------------------------------
naginterfaces.base.opt.handle_solve_socp_ipm: Final primal objective value -8.505797E-02
naginterfaces.base.opt.handle_solve_socp_ipm: Final dual objective value -8.505797E-02
naginterfaces.base.opt.handle_solve_socp_ipm: Absolute primal infeasibility 2.869531E-09
naginterfaces.base.opt.handle_solve_socp_ipm: Relative primal infeasibility 3.186829E-10
naginterfaces.base.opt.handle_solve_socp_ipm: Absolute dual infeasibility 1.765929E-09
naginterfaces.base.opt.handle_solve_socp_ipm: Relative dual infeasibility 1.471607E-10
naginterfaces.base.opt.handle_solve_socp_ipm: Absolute complementarity gap 1.266105E-09
naginterfaces.base.opt.handle_solve_socp_ipm: Relative complementarity gap 8.044616E-10
naginterfaces.base.opt.handle_solve_socp_ipm: Iterations 11これで最適なポートフォリオとそれに対応するリターンを表示できます。
# 最適ポートフォリオ
slt.x[0:n_assets]array([1.93322779e-10, 3.97165177e-01, 3.35524516e-01, 7.15459782e-10,
1.16762817e-01, 1.15110232e-10, 8.27684160e-02, 6.77790766e-02])# 最適期待収益率
r.dot(slt.x[0:n_assets])0.08505797359531912結論
このノートブックでは、nAGライブラリを使用して、SOCPを介した確率制約付きポートフォリオ最適化をモデル化し解決する方法を示しました。SOCPがポートフォリオ最適化において広く使用されているのは、その柔軟性と汎用性により、上記の確率制約だけでなく、レバレッジ制約、回転率制約、最大ポジション制約、トラッキングエラー制約など、さまざまな種類の制約を持つ多様な問題を扱えるためです。詳細については、を参照してください。
参考文献
[1] Alizadeh Farid and Goldfarb Donald, ``Second-order cone programming’’, Mathematical programming, vol. 95, number 1, pp. 3–51, 2003.
[2] Lobo Miguel Sousa, Vandenberghe Lieven, Boyd Stephen et al., ``Applications of second-order cone programming’’, Linear algebra and its applications, vol. 284, number 1-3, pp. 193–228, 1998.
[3] Numerical Algorithms Group, ``nAG documentation’’, 2019. online