混合整数二次計画法ポートフォリオ最適化問題: 問題ベース

ライブスクリプトを開く

この例では、問題ベースアプローチを使用して、混合整数二次計画法 (MIQP) のポートフォリオ最適化問題の解法を示します。考え方としては、MIQP 問題を局所的に近似する一連の混合整数線形計画法 (MILP) の解を反復的に求めるというものです。ソルバーベースのアプローチについては、混合整数二次計画法ポートフォリオ最適化問題: ソルバーベースを参照してください。

問題の概要

Markowitz が示したように、多くのポートフォリオ最適化問題は二次計画問題として表されます ("Portfolio Selection," J. Finance Volume 7, Issue 1, pp. 77-91, March 1952)。N 件の項目から成る一連の資産があり、ポートフォリオを選ぶとします。ここで $x (i)$ は資産 $i$ における投資割合です。各資産の平均収益のベクトル $r$ および収益の共分散行列 $Q$ がわかっている場合、与えられた危険回避のレベル $λ$ に対し、リスク調整後の期待収益を最大化します。

$\max_{x} (r^{T} x - λ x^{T} Q x) .$

quadprog ソルバーを使用してこの二次計画問題を解きます。しかし、そのままの二次計画問題に加えて、ポートフォリオに次のようなさまざまな制限をかけることが必要な場合があります。

ポートフォリオに含まれる資産が M 個を超えないものとする。ここで M <= N。
ポートフォリオが少なくとも m 個の資産を含むとする。ここで 0 < m <= M。
"半連続的な" 制約があるものとする。これは、 $x (i) = 0$ であるか、 $f m i n > 0$ および $f m a x \geq f m i n$ のいくつかの固定割合について $f m i n \leq x (i) \leq f m a x$ であることを意味します。

quadprog にこれらの制約を含めることはできません。難しいのは、制約の離散的な性質です。さらに、混合整数線形計画法ソルバーは離散的な制約は処理できますが、二次目的関数は扱いません。

この例では、制約を満たすと共に、しだいに二次の目的関数を近似するようになる一連の MILP 問題を作成します。この手法はこの例については機能しますが、他の問題や制約のタイプには当てはまらない場合があります。

まず制約をモデル化します。

離散的な制約のモデル化

$x$ は資産配分の割合を表すベクトルであり、各 $i$ について $0 \leq x (i) \leq 1$ です。ポートフォリオに含まれる資産の個数をモデル化するため、 $x (i) = 0$ の場合は $v (i) = 0$ 、 $x (i) > 0$ の場合は $v (i) = 1$ になるように指標変数 $v$ を設定します。この制限を満たす変数を得るため、ベクトル $v$ をバイナリ変数に設定し、線形制約を課します。

$v (i) f m i n \leq x (i) \leq v (i) f m a x .$

これらの不等式によって、 $x (i)$ と $v (i)$ はまったく同時にゼロになります。また、 $x (i) > 0$ の場合は常に、 $f m i n \leq x (i) \leq f m a x$ になります。

また、ポートフォリオ内の資産数への制約を実現させるため、線形制約を課します。

$m \leq \sum_{i} v (i) \leq M .$

目的関数および連続的線形近似

最初に定式化したように、目的関数を最大化しようと試みます。ところが、Optimization Toolbox™ のソルバーではすべて最小化されます。そこで、目的関数の負の値を最小化するように問題を定式化します。

$\min_{x} λ x^{T} Q x - r^{T} x .$

この目的関数は非線形です。MILP ソルバーには線形目的関数が必要です。この問題を線形目的関数と非線形な制約をもつ問題に再度定式化する標準の手法があります。二次の項を表すためにスラック変数 $z$ を導入します。

$\min_{x, z} λ z - r^{T} x such that x^{T} Q x - z \leq 0, z \geq 0 .$

MILP 近似を反復的に解く際、新規の線形制約を組み込み、各制約によって、現在の点の近くにある非線形な制約を局所的に近似します。具体的には、 $x_{0}$ が定数ベクトル、 $δ$ が変数ベクトルである $x = x_{0} + δ$ について、制約の 1 次テイラー近似は次のようになります。

$x^{T} Q x - z = x_{0}^{T} Q x_{0} + 2 x_{0}^{T} Q δ - z + O (| δ |^{2}) .$

$δ$ を $x - x_{0}$ で置き換えると、次のようになります。

$x^{T} Q x - z = - x_{0}^{T} Q x_{0} + 2 x_{0}^{T} Q x - z + O (| x - x_{0} |^{2}) .$

各中間解 $x_{k}$ について、上式の線形部分として $x$ と $z$ に新しい線形制約を導入します。

$- x_{k}^{T} Q x_{k} + 2 x_{k}^{T} Q x - z \leq 0 .$

これは $A x \leq b$ という形式になります。ここで、 $A = 2 x_{k}^{T} Q$ 、 $z$ 項の乗数は $- 1$ 、 $b = x_{k}^{T} Q x_{k}$ です。

問題に新規の線形制約を追加するこの方法は、切除平面法と呼ばれます。詳細については、J. E. Kelley, Jr. "The Cutting-Plane Method for Solving Convex Programs." J. Soc.Indust.Appl.Math.Vol. 8, No. 4, pp. 703-712, December, 1960 を参照してください。

MATLAB® での問題の定式化

最適化問題を表すには、次を行います。

変数が表す内容の決定
これらの変数の下限と上限の明示
線形の等式および不等式の指定

問題のデータを読み込みます。このデータには、ベクトル r に含まれる 225 個の期待される収益と、225 行 225 列の行列 Q に含まれる収益の共分散があります。データは、「ポートフォリオ最適化問題に対する二次計画法の使用」の例と同じです。

load port5
r = mean_return;
Q = Correlation .* (stdDev_return * stdDev_return');

資産の数を N に設定します。

N = length(r);

問題の変数、制約および目的の作成

資産配分の割合を表す連続変数 xvars、関連付けられた xvars がゼロと厳密に正のどちらであるかを示すバイナリ変数 vvars、および正のスカラーである変数 $z$ を表す zvar を作成します。

xvars = optimvar('xvars',N,1,'LowerBound',0,'UpperBound',1);
vvars = optimvar('vvars',N,1,'Type','integer','LowerBound',0,'UpperBound',1);
zvar = optimvar('zvar',1,'LowerBound',0);

問題において 2N+1 個のすべての変数の下限はゼロです。変数 xvars および yvars の上限は 1 で、zvar には上限はありません。

解における資産の数を 100 ～ 150 に設定します。この制約を次の形式で問題に組み込みます。

$m \leq \sum_{i} v (i) \leq M,$

2 つの線形制約で記述するとこうなります。

$\sum_{i} v (i) \leq M$

$\sum_{i} v (i) \geq m .$

M = 150;
m = 100;
qpprob = optimproblem('ObjectiveSense','maximize');
qpprob.Constraints.mconstr = sum(vvars) <= M;
qpprob.Constraints.mconstr2 = sum(vvars) >= m;

半連続的な制約を含めます。それぞれの資産タイプについて、非ゼロの最小の資産の割合を 0.001 に、また、最大の割合を 0.05 にします。

fmin = 0.001;
fmax = 0.05;

不等式 $x (i) \leq f m a x (i) * v (i)$ および $f m i n (i) * v (i) \leq x (i)$ を含めます。

qpprob.Constraints.fmaxconstr = xvars <= fmax*vvars;
qpprob.Constraints.fminconstr = fmin*vvars <= xvars;

ポートフォリオは 100% 運用であるという制約を含めます。これは、 $\sum x_{i} = 1$ を意味します。

qpprob.Constraints.allin = sum(xvars) == 1;

危険回避係数 $λ$ を 100 に設定します。

lambda = 100;

目的関数 $r^{T} x - λ z$ を定義し、問題に含めます。

qpprob.Objective = r'*xvars - lambda*zvar;

問題を解く

問題を反復的に解くために、現在の制約をもつ問題を解くことから始めます。現在の制約については、いかなる線形化も実施していません。

options = optimoptions(@intlinprog,'Display','off'); % Suppress iterative display
[xLinInt,fval,exitFlagInt,output] = solve(qpprob,'options',options);

反復の停止条件を決めます。スラック変数 $z$ が真の二次値の 0.01% 以内になったときに反復を止めます。

thediff = 1e-4;
iter = 1; % iteration counter
assets = xLinInt.xvars;
truequadratic = assets'*Q*assets;
zslack = xLinInt.zvar;

プロットするために二次の真値とスラック変数の計算履歴を保存します。反復が正しい解に収束するように、許容誤差を既定より小さくします。

history = [truequadratic,zslack];

options = optimoptions(options,'LPOptimalityTolerance',1e-10,'RelativeGapTolerance',1e-8,...
                      'ConstraintTolerance',1e-9,'IntegerTolerance',1e-6);

二次の値およびスラック変数の値を計算します。それらが異なる場合は、別の線形制約を追加して再度解を求めます。

新しい各線形制約 $A x \leq b$ は次の線形近似から得られます。

$- x_{k}^{T} Q x_{k} + 2 x_{k}^{T} Q x - z \leq 0 .$

新規の解が求まったら、元の解と新規解の中間にある線形制約を使用します。線形制約を含んでいるこのヒューリスティックな方法のほうが、単純に新規解を用いるよりも高速になる場合があります。この中間のヒューリスティックな方法の代わりに解を使う場合は、下記の "Midway" の行にコメントし、その次の行のコメントを解除します。

while abs((zslack - truequadratic)/truequadratic) > thediff % relative error
    constr = 2*assets'*Q*xvars - zvar <= assets'*Q*assets;
    newname = ['iteration',num2str(iter)];
    qpprob.Constraints.(newname) = constr;
    % Solve the problem with the new constraints
    [xLinInt,fval,exitFlagInt,output] = solve(qpprob,'options',options);
    assets = (assets+xLinInt.xvars)/2; % Midway from the previous to the current
%     assets = xLinInt(xvars); % Use the previous line or this one
    truequadratic = xLinInt.xvars'*Q*xLinInt.xvars;
    zslack = xLinInt.zvar;
    history = [history;truequadratic,zslack];
    iter = iter + 1;
end

解と収束率の検証

スラック変数と目的関数の二次の部分の計算履歴をプロットしてどのように収束するかを確認します。

plot(history)
legend('Quadratic','Slack')
xlabel('Iteration number')
title('Quadratic and linear approximation (slack)')

MILP の解はどのような品質でしょうか。output 構造体に、その情報が含まれます。解における、内部で計算された目的関数の限界値間のギャップの絶対値を検証します。

disp(output.absolutegap)

ギャップの絶対値はゼロであり、MILP の解は正確であることを示しています。

最適な資産配分をプロットします。中間位置での更新を使用すると、assets は制約を満たさない可能性があるため、assets ではなく、xLinInt.xvars を使用します。

bar(xLinInt.xvars)
grid on
xlabel('Asset index')
ylabel('Proportion of investment')
title('Optimal asset allocation')

非ゼロの資産配分のすべてが半連続的な範囲 $f m i n = 0.001$ と $f m a x = 0.05$ の間にあることが簡単に確認できます。

非ゼロの資産がいくつあるでしょうか。制約は、非ゼロの資産の数が 100 ～ 150 です。

sum(xLinInt.vvars)

ans = 100

この資産配分により期待される収益、そしてリスク調整された収益値はどれほどでしょうか。

fprintf('The expected return is %g, and the risk-adjusted return is %g.\n',...
    r'*xLinInt.xvars,fval)

The expected return is 0.000595107, and the risk-adjusted return is -0.0360382.

特にポートフォリオ最適化用に設計された機能を Financial Toolbox® で使用すると、より詳細な分析が可能です。ポートフォリオクラスを使用して半連続および濃度の制約を直接処理する方法を示す例については、Portfolio Optimization with Semicontinuous and Cardinality Constraints (Financial Toolbox)を参照してください。