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

問題の概要

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

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

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

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

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

離散的な制約のモデル化

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

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

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

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

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

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

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

$$\delta$$ を $$x-x_0$$ で置き換えると、次のようになります。

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

これは $$Ax\le b$$ という形式になります。ここで、$$A=2x_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® での問題の定式化

intlinprog ソルバー向けに問題を表すためには以下を行う必要があります。

最初の $$N$$ 個の変数は $$x$$ ベクトルを表し、次の $$N$$ 個の変数は 2 値の $$v$$ ベクトルを表し、最後の変数はスラック変数 $$z$$ を表します。問題には $$2N+1$$ 個の変数が含まれています。

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

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

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

N = length(r);

変数のインデックスを設定します

xvars = 1:N;
vvars = N+1:2*N;
zvar = 2*N+1;

問題において 2N+1 個のすべての変数の下限はゼロです。最初の 2N 個の変数の上限は 1 で、最後の変数には上限はありません。

lb = zeros(2*N+1,1);
ub = ones(2*N+1,1);
ub(zvar) = Inf;

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

形式 $$Ax\le b$$ の 2 つの線形制約を記述すると、次のようになります。

M = 150;
m = 100;
A = zeros(1,2*N+1); % Allocate A matrix
A(vvars) = 1; % A*x represents the sum of the v(i)
A = [A;-A];
b = zeros(2,1); % Allocate b vector
b(1) = M;
b(2) = -m;

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

fmin = 0.001;
fmax = 0.05;

線形不等式として、不等式 $$x(i)\le fmax(i)*v(i)$$ および $$fmin(i)*v(i)\le x(i)$$ を含めます。

Atemp = eye(N);
Amax = horzcat(Atemp,-Atemp*fmax,zeros(N,1));
A = [A;Amax];
b = [b;zeros(N,1)];
Amin = horzcat(-Atemp,Atemp*fmin,zeros(N,1));
A = [A;Amin];
b = [b;zeros(N,1)];

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

Aeq = zeros(1,2*N+1); % Allocate Aeq matrix
Aeq(xvars) = 1;
beq = 1;

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

lambda = 100;

目的関数 $$\lambda z-r^{T}x$$ をベクトルとして定義します。変数 $$v$$ の乗数にはゼロを含めます。

f = [-r;zeros(N,1);lambda];

問題を解く

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

options = optimoptions(@intlinprog,'Display','off'); % Suppress iterative display
[xLinInt,fval,exitFlagInt,output] = intlinprog(f,vvars,A,b,Aeq,beq,lb,ub,options);

反復の停止条件を決めます。スラック変数 $$z$$ が真の二次値の 0.01% 以内になったときに反復を止めます。制約が累積されても問題が厳密に実行可能であることを確保するために、制約の許容誤差を既定より小さく設定します。

thediff = 1e-4;
iter = 1; % iteration counter
assets = xLinInt(xvars); % the x variables
truequadratic = assets'*Q*assets;
zslack = xLinInt(zvar); % slack variable value
options.ConstraintTolerance = 1e-9;

プロットするために二次の真値とスラック変数の計算履歴を保存します。

history = [truequadratic,zslack];

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

ツールボックスの構文にある、それぞれの新規の線形制約 $$Ax\le b$$ は次の線形近似から得られます。

$$A=2x_k^{T}Q$$ の新しい行、$$b=x_k^{T}Q{x}_k$$ の新しい要素、および $$A$$ の -1 係数で表される $$z$$ 項があるのがわかります。

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

while abs((zslack - truequadratic)/truequadratic) > thediff % relative error
    newArow = horzcat(2*assets'*Q,zeros(1,N),-1); % Linearized constraint
    rhs = assets'*Q*assets;                       % right hand side of the linearized constraint
    A = [A;newArow];
    b = [b;rhs];
    % Solve the problem with the new constraints
    [xLinInt,fval,exitFlagInt,output] = intlinprog(f,vvars,A,b,Aeq,beq,lb,ub,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)

     0

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

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

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

非ゼロの資産配分のすべてが半連続的な範囲 $$fmin=0.001$$ と $$fmax=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)を参照してください。