nAGライブラリを使用したMILPによるポートフォリオ最適化

混合整数線形計画法（MILP）モデルは、一部または全ての変数が整数に制約される線形計画法モデルの拡張です。整数変数を扱う能力により、これは非常に強力なツールとなり、膨大な数の産業に応用されています。ここでは、Benati and Rizzi (2007)による最適な平均/バリューアットリスク（VaR）ポートフォリオ最適化問題のMILPモデルを考察します。これは、リターンを最大化しながらリスクを最小化することを目的とする古典的なマーコビッツモデルの拡張です。ただし、この場合、分散リスク指標がVaRに置き換えられています。

VaRは、ポートフォリオ管理で広く使用されているリスク指標です：特定の信頼水準で定義された期間内の最大潜在損失を定量化します。具体的には、VaRは単にリターン分布関数の\(\alpha\)分位数です。標準偏差や期待ショートフォールなどの他のリスク指標とは異なり、VaRは下方リスクの明確で直感的な評価を提供し、重大な損失を経験する確率に関心のある投資家にとって好ましい選択肢となっています。VaRの数学的特性はやや魅力的ではありませんが（区分線形関数であり凸ではありません）、その解釈の簡潔さにより、VaRは依然として人気のあるツールです。

この場合、VaRはノンパラメトリック法を用いて計算されます - これは、ポートフォリオリターンの分布に関する仮定を行わないことを意味します。歴史的データにパラメトリックモデル（例：正規分布やt分布）を当てはめる代わりに、歴史的シミュレーションを使用します。この方法は、歴史的リターンデータを最悪から最良の順に並べ替えます。そして、\(\alpha\)分位数は、左側に\(\alpha\)パーセントのデータを持つ観測値の位置によって推定されます。

このノートブックでは、この問題の「最小リスク/固定リターン」実装を解きます。この問題には、VaR制約、カーディナリティ制約、半連続制約、最小リターン制約、予算（全額投資）制約、およびロングオンリー制約が含まれています。

# モジュールをインポート
from naginterfaces.base import utils
from naginterfaces.library import opt, mip
import numpy as np
import time

このモデルでは、いくつかのパラメータを設定する必要があります。

過去の観測値 \(i \in I=\{1,\dots,T\}\) と資産 \(j \in J=\{1,\dots,K\}\) があります。

\(r^*\) を受け入れ可能な最小期待リターンとし、\(r^{\textrm{Min}}\) は市場で観測可能な最小リターンとします。また、\(r^{\textrm{VaR}}\) を選択します：これは意思決定者がリスクを制御するために設定するパラメータです。リターンが \(r^{\textrm{VaR}}\) 未満となる確率が \(\alpha^{\textrm{VaR}}\) 以下であるポートフォリオのみを受け入れます。さらに、各観測値 \(x_i\) に関連する確率 \(p_i\) があります。これらの確率は、過去の実現値 \(i\) の発生を表します。説明のために、各観測値に等しい確率を設定し、リターンデータを人工的に生成します。

# パラメータを設定
n_assets = 300 # K
n_periods = 20 # T
r_star = 0.05
r_min = -1
r_var = 0.05
prob = 1/n_periods

# 各資産の期待リターンの合成データ生成
np.random.seed(0)
expected_returns = 0.25 * np.random.randn(n_periods, n_assets) 

次に、モデルで使用される変数を定義します： * 変数 \(\lambda_j\) は資産 \(j\) に配分される富の割合です。 * 変数 \(x_i\) は時間 \(i\) におけるポートフォリオの観測リターンを表します。 * 変数 \(y_i\) はバイナリ（\(x_i\) に関連）で、VaR制約のモデリングに使用されます。 * 変数 \(z_j\) はバイナリ（\(\lambda_j\) に関連）で、カーディナリティと半連続制約のモデリングに使用されます。 * \(\alpha^{\textrm{VaR}}\) はVaRに関連する確率です。

# 以下の変数があります：asset_weights[n_assets] + binary_y[n_periods] + observed_return[n_periods] + binary_z[n_assets] + aVaR[1]
n_vars = n_assets + n_periods + n_periods + n_assets + 1

# 各変数セットのインデックス集合を作成します:
idx_asset_w = np.arange(1, n_assets+1, dtype=int)
idx_bin_y = np.arange(n_assets+1, n_assets+n_periods+1, dtype=int)
idx_obs_ret = np.arange(n_assets+n_periods+1, n_assets+2*n_periods+1, dtype=int)
idx_bin_z = np.arange(n_assets+2*n_periods+1, 2*n_assets+2*n_periods+1, dtype=int)
idx_avar = [n_vars]

# 問題ハンドルを初期化します:
handle = opt.handle_init(nvar=n_vars)

# 2値変数を設定:
opt.handle_set_property(handle=handle, ptype='Bin', idx=idx_bin_y)
opt.handle_set_property(handle=handle, ptype='Bin', idx=idx_bin_z)

制約条件:

これはMILPモデルであるため、すべての線形制約を handle_set_linconstr() を使用して設定します。

実装する最初の制約は、最適なポートフォリオが最小許容ポートフォリオ期待収益率よりも大きくなければならないことを強制します: \[\sum_{i=1}^{T} p_i x_i \geq r^*.\]

ones_periods = np.ones(n_periods, dtype=int)
opt.handle_set_linconstr(
    handle=handle,
    bl=r_star,
    bu=1.e20,
    irowb=ones_periods,
    icolb=idx_obs_ret,
    b=[prob]*n_periods
);

次に、各時期において\(x_i\)が各資産に配分された富の割合とその資産の期待収益率の結果であることを強制する必要があります：

\[\sum_{j=1}^{K} \lambda_j r_{ij} = x_i \quad \forall \, i=1,\dots,T.\]

これをモデルに入力するには、若干の並べ替えが必要です - 変数\(\lambda\)と\(x\)を積み重ね、並べ替えを行列表記で表現します：

\[\begin{equation*} \begin{bmatrix} r_{i1} & \dots & r_{iK} & -1 \end{bmatrix} \begin{bmatrix} \lambda_1 \\ \vdots \\ \lambda_K \\ x_i \end{bmatrix} = 0 \quad \forall \, i = 1,\dots,T. \end{equation*}\]

ones_assets_1 = np.ones(n_assets+1, dtype=int)
for i in range(n_periods):
    i_obs_ret = n_assets + n_periods + 1 + i
    returns = expected_returns[i, :]
    returns_x = np.append(returns, -1.)
    idx_col = np.append(idx_asset_w, i_obs_ret)
    opt.handle_set_linconstr(
        handle=handle,
        bl=0.,
        bu=0.,
        irowb=ones_assets_1,
        icolb=idx_col,
        b=returns_x
    )

次に、VaRが閾値を下回るポートフォリオの選択を防ぐ制約を追加したいと思います。これを行うために、各\(x_i\)に関連するバイナリ変数\(y_i\)を導入します。そうすると、VaR制約は以下のようにモデル化できます：

\[r^{\textrm{Min}} + (r^{\textrm{VaR}} - r^{\textrm{Min}})y_i \leq x_i \quad \forall \, i=1,\dots,T,\] \[\sum_{i=1}^{T} p_i(1-y_i) \leq \alpha^{\textrm{VaR}}.\]

最初の制約は、\(x_i\)が\(r^{\textrm{VaR}}\)未満の場合、\(y_i\)が\(0\)に等しくなることを強制します。2番目の制約では、これは\(1 - y_i = 1 - 0 = 1\)に対応し、VaR閾値未満のリターンを持つ期間\(i\)の確率の合計につながります。この確率が\(\alpha^{\textrm{VaR}}\)より大きい場合、実行不可能なポートフォリオとなります。

行列表記に並べ替えると、最初の制約は以下のようになります：

\[\begin{equation*} \begin{bmatrix} (r^{\textrm{Min}}-r^{\textrm{VaR}}) & 1 \end{bmatrix} \begin{bmatrix} y_i \\ x_i \end{bmatrix} \geq r^{\textrm{Min}} \quad \forall \, i = 1,\dots,T, \end{equation*}\]

そして2番目の制約は以下のようになります：

\[\begin{equation*} \begin{bmatrix} p_{1} & \dots & p_{T} & 1 \end{bmatrix} \begin{bmatrix} y_1 \\ \vdots \\ y_T \\ \alpha^{\textrm{VaR}} \end{bmatrix} \geq 1. \end{equation*}\]

y_coef = r_min - r_var
for i in range(n_periods):
    y_i = n_assets + 1 + i
    x_i = n_assets + n_periods + 1 + i     
    opt.handle_set_linconstr(
        handle=handle,
        bl=r_min,
        bu=1.e20,
        irowb=[1, 1],
        icolb=[y_i, x_i],
        b=[y_coef, 1.]
    )

ones_periods_1 = np.ones(n_periods+1, dtype=int)
idx_bin_avar = np.append(idx_bin_y, idx_avar)
opt.handle_set_linconstr(
    handle=handle,
    bl=1.0,
    bu=1.e20,
    irowb=ones_periods_1,
    icolb=idx_bin_avar,
    b=np.append([prob]*n_periods, 1.)
);

全額投資制約:

\[\sum_{j=1}^{K} \lambda_j = 1.\]

ones_assets_float = np.ones(n_assets, dtype=float)
ones_assets = np.ones(n_assets, dtype=int)
opt.handle_set_linconstr(
    handle=handle,
    bl=1.0,
    bu=1.0,
    irowb=ones_assets,
    icolb=idx_asset_w,
    b=ones_assets_float
);

次に、資産の最適配分を5%から70%の間に制限する半連続制約を追加したいと思います: \[\lambda_j \in 0 \cup [0.05, 0.7] \quad \forall \, j = 1,\dots, K.\]

バイナリ変数 \(z_j\) を使用すると、半連続制約は以下のように表現できます:

\[0.05 \cdot z_j \leq \lambda_j \leq 0.7 \cdot z_j \quad \forall \, j=1,\dots, K.\]

これらの制約は、\(z_j\) が \(1\) に等しい場合、\(\lambda_j\) は \(0.05\) から \(0.7\) の間の任意の値を取ることができることを強制します。\(z_j\) が \(0\) に等しい場合、\(\lambda_j\) も \(0\) に等しくなければなりません。

for j in range(n_assets):
    lambda_j = 1 + j
    z_j = n_assets + 2*n_periods + 1 + j
    opt.handle_set_linconstr(
        handle=handle,
        bl=-1.e20,
        bu=0.0,
        irowb=[1, 1],
        icolb=[lambda_j, z_j],
        b=[1., -0.7]
    )
   
    opt.handle_set_linconstr(
        handle=handle,
        bl=0.0,
        bu=1.e20,
        irowb=[1, 1],
        icolb=[lambda_j, z_j],
        b=[1., -0.05]
    )

次に、最適ポートフォリオに資産の10%を含めたいので、基数制約を追加します： \[||\lambda||_0 \leq 0.1 \cdot K.\]

これは、各\(\lambda_j\)に対してバイナリ変数\(z_j\)を追加することでモデル化できます。\(z_j\)が0の場合、資産\(\lambda_j\)が割り当てられていないことを示し、\(z_j\)が1の場合、資産\(\lambda_j\)が割り当てられたことを示します。そして、基数制約はバイナリ変数の合計によってモデル化できます：

\[\sum_{j=1}^K z_j \leq 0.1 \cdot K.\]

opt.handle_set_linconstr(
    handle=handle,
    bl=-1.e20,
    bu=0.1*n_assets,
    irowb=ones_assets,
    icolb=idx_bin_z,
    b=ones_assets_float
);

空売りを防ぐための境界制約：

\[\lambda_j \geq 0 \quad \forall \, j=1,\dots,K.\]

# 境界制約を設定:
zeros_asset_w = np.zeros(n_assets)
neg_lrg_bnd = np.full((n_periods + n_periods + n_assets + 1), -1.e20, dtype=float).ravel()
bl = np.hstack([zeros_asset_w, neg_lrg_bnd])
opt.handle_set_simplebounds(
    handle,
    bl=bl,
    bu=[1.e20]*n_vars,
)

目的関数:

リスクを最小化したいと考えています:

\[\begin{equation*} \min_{\alpha^{\textrm{VaR}},\lambda,x,y} \alpha^{\textrm{VaR}} \end{equation*}\]

目的関数は線形ですが、疎であるためhandle_set_quadobj()を使用します。

# 目的関数を設定する:
opt.handle_set_quadobj(
        handle,
        idxc=idx_avar,
        c=[1.0],
    )

最適平均/バリューアットリスクモデル

ここで、完全なモデルを指定し、MILP ソルバーを使用して解くことができます。

\[\begin{equation*} \begin{aligned} &\min_{\alpha^{\textrm{VaR}},\lambda,x,y,z} &&\alpha^{\textrm{VaR}} \\ &\textrm{subject to } &&\sum_{i=1}^{T} p_i x_i \geq r^*, \\ & &&\sum_{j=1}^{K} \lambda_j r_{ij} = x_i \quad \forall \, i=1,\dots,T, \\ & &&r^{\textrm{Min}} + (r^{\textrm{VaR}} - r^{\textrm{Min}})y_i \leq x_i \quad \forall \, i=1,\dots,T, \\ & &&\sum_{i=1}^{T} p_i(1-y_i) \leq \alpha^{\textrm{VaR}}, \\ & &&\sum_{j=1}^{K} \lambda_j = 1, \\ & &&\sum_{j=1}^K z_j \leq 0.1 \cdot K, \\ & && 0.05 \cdot z_j \leq \lambda_j \leq 0.7 \cdot z_j \quad \forall \, j = 1,\dots, K, \\ & &&y_i \in \{0,1\} \quad \forall \, i=1,\dots,T, \\ & &&z_j \in \{0,1\} \quad \forall \, j=1,\dots,K, \\ & &&\lambda_j \geq 0 \quad \forall \, j=1,\dots,K. \end{aligned} \end{equation*}\]

# アルゴリズムのオプションをいくつか設定します:
for option in [
        'Print Level = 1',
]:
    opt.handle_opt_set(handle, option)

# 省略された反復出力のために明示的なI/Oマネージャーを使用する：
iom = utils.FileObjManager(locus_in_output=False)

# タイマーを開始:
start_solve = time.time()

# MILPソルバーを呼び出す:
mip.handle_solve_milp(handle, io_manager=iom)

# 計算時間を表示:
end = time.time()
print(f"\n Computation time: {end-start_solve:.3f} seconds")

# ハンドルを破棄する:
opt.handle_free(handle)

 H02BK, Solver for MILP problems
 Begin of Options
     Print File                    =                   9     * d
     Print Level                   =                   1     * U
     Print Options                 =                 Yes     * d
     Print Solution                =                  No     * d
     Monitoring File               =                  -1     * d
     Monitoring Level              =                   4     * d

     Infinite Bound Size           =         1.00000E+20     * d
     Task                          =            Minimize     * d
     Time Limit                    =         1.00000E+06     * d

     Milp Presolve                 =                 Yes     * d
     Milp Random Seed              =                   0     * d
     Milp Feasibility Tol          =         1.00000E-06     * d
     Milp Rel Gap                  =         1.00000E-04     * d
     Milp Abs Gap                  =         1.00000E-06     * d
     Milp Small Matrix Value       =         1.00000E-09     * d
     Milp Detect Symmetry          =                 Yes     * d
     Milp Max Nodes                =          2147483647     * d
 End of Options

 Status: converged, an optimal solution found
 Final primal objective value -5.551115E-17
 Final dual objective bound    0.000000E+00

 Computation time: 0.405 seconds

この例では、641個の変数（そのうち320個がバイナリ）と644個の制約を持つモデルを解きました。これは0.41秒で効率的に解かれました。

MILPを使用することで、複雑な制約のモデリングが可能になり、問題に柔軟性と洗練さを加えることができます。このポートフォリオ最適化の例で使用されたMILPソルバーは、nAGライブラリ最適化モデリングスイートで利用可能です。詳細はこちらをご覧ください。

参考文献

Benati, S., & Rizzi, R. (2007). A mixed integer linear programming formulation of the optimal mean/value-at-risk portfolio problem. European Journal of Operational Research, 176(1), 423-434.