exponenta event banner

Портфолио квадратичного программирования со смешанным целым: на основе решателя

В этом примере показано, как решить проблему оптимизации портфеля Mike-Integer Quadratic Programming (MIQP) с помощью intlinprog Решатель смешанного целочисленного линейного программирования (MILP). Идея состоит в итеративном решении последовательности задач MILP, которые локально аппроксимируют задачу MIQP. Для получения информации о подходе, основанном на проблемах, см. раздел Оптимизация портфеля квадратичного программирования со смешанными целыми значениями: на основе проблем.

Структура проблемы

Как показал Марковиц («Portfolio Selection», J. Finance Volume 7, Issue 1, pp. 77-91, March 1952), многие задачи оптимизации портфеля можно выразить как задачи квадратичного программирования. Предположим, что у вас есть набор N Если вы знаете вектор r средней доходности каждого актива и ковариационную матрицу Q доходности, то для данного уровня отклонения от риска λ вы максимизируете скорректированную на риск ожидаемую доходность:

maxx (rTx-λ xTQx).

quadprog решатель решает эту проблему квадратичного программирования. Однако в дополнение к проблеме простого квадратичного программирования можно ограничить портфолио различными способами, такими как:

  • Имея не более M активы в портфеле, где M <= N.

  • Имея по крайней мере m активы в портфеле, где 0 < m <= M.

  • Имеет полунепрерывные ограничения, означающие либо x (i) = 0, либо fmin≤x (i) ≤fmax для некоторых фиксированных дробей fmin > 0 и fmax≥fmin.

Нельзя включать эти ограничения в quadprog. Трудность заключается в дискретном характере ограничений. Кроме того, в то время как решатель линейного программирования со смешанным целым числом intlinprog обрабатывает дискретные ограничения, не обращается к квадратичным объективным функциям.

Этот пример строит последовательность задач MILP, которые удовлетворяют ограничениям, и которые все больше аппроксимируют квадратичную целевую функцию. Хотя этот метод работает в этом примере, он может не применяться к другим типам проблем или ограничений.

Начните с моделирования зависимостей.

Моделирование дискретных ограничений

x - вектор дробей распределения основных средств, с 0≤x (i) ≤1 для каждого из них. Для моделирования количества основных средств в портфеле необходимы переменные индикатора v, так что v (i) = 0, когда x (i) = 0, и v (i) = 1, когда x (i) > 0. Чтобы получить переменные, удовлетворяющие этому ограничению, задайте вектор v как двоичную переменную и наложите линейные ограничения

v (i) fmin≤x (i) ≤v (i) fmax.

Эти неравенства обеспечивают, что x (i) и v (i) равны нулю точно в одно и то же время, и они также обеспечивают, что fmin≤x (i) ≤fmax всякий раз, когда x (i) > 0.

Кроме того, для обеспечения соблюдения ограничений на количество активов в портфеле, наложить линейные ограничения

m≤∑iv (i) ≤M.

Целевые и последовательные линейные приближения

В соответствии с первой формулировкой, вы пытаетесь максимизировать целевую функцию. Однако все решатели Optimization Toolbox™ минимизируются. Так что сформулируйте проблему как минимизирующую негатив цели:

minxλ xTQx-rTx.

Эта целевая функция нелинейна. intlinprog Решателю MILP требуется линейная целевая функция. Существует стандартный метод переформулировать эту проблему в задачу с линейными объективными и нелинейными ограничениями. Введите переменную провисания z для представления квадратичного члена.

minx,  zλ z-rTx такие , что xTQx-z≤0,z≥0.

При итеративном решении аппроксимаций MILP включаются новые линейные ограничения, каждое из которых аппроксимирует нелинейное ограничение локально вблизи текущей точки. В частности, для x = x0 + δ, где x0 - постоянный вектор, а δ - переменный вектор, аппроксимация Тейлора первого порядка к ограничению равна

xTQx-z = x0TQx0 + 2x0TQδ-z + O (| δ | 2).

Замена δ на x-x0 дает

xTQx-z = -x0TQx0 + 2x0TQx-z + O (| x-x0 | 2).

Для каждого промежуточного решения xk вводится новое линейное ограничение в x и z в качестве линейной части выражения выше:

- xkTQxk+2xkTQx-z≤0.

Имеет вид Ax≤b, где A = 2xkTQ, имеется множитель -1 для z-члена и b = xkTQxk.

Этот метод добавления новых линейных ограничений к задаче называется методом секущей плоскости. Подробнее см. J. E. Kelley, Jr. «Метод секущей плоскости для решения выпуклых программ». Дж. Соц. Индуст. прим. мат. т. 8, № 4, стр. 703-712, декабрь 1960.

Формулировка проблемы MATLAB ®

Чтобы выразить проблемы для intlinprog решатель, необходимо выполнить следующие действия.

  • Решите, что представляют ваши переменные

  • Выражать нижние и верхние границы с точки зрения этих переменных

  • Дать матрицы линейного равенства и неравенства

Первые N переменных представляют вектор x, следующие N переменных представляют двоичный v-вектор, а окончательная переменная представляет переменную z slack. В проблеме есть переменные 2N + 1.

Загрузите данные для проблемы. Эти данные имеют 225 ожидаемых результатов в векторе r и ковариация возвратов в матрице 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 переменные равны единице, а последняя переменная не имеет верхней границы.

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

Установите количество основных средств в решении в диапазоне от 100 до 150. Включить это ограничение в проблему в форме, а именно

m≤∑iv (i) ≤M,

путем записи двух линейных ограничений вида Ax≤b:

∑iv (i) ≤M

∑i-v (i) ≤-m.

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) ≤fmax (i) * v (i) и fmin (i) * v (i) ≤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%, что означает ∑xi=1.

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

Установите коэффициент отклонения от риска λ в значение 100.

lambda = 100;

Определите целевую функцию λ z-rTx как вектор. Включите нули для множителей переменных 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 = optimoptions(options,'LPOptimalityTolerance',1e-10,'RelativeGapTolerance',1e-8,...
                      'ConstraintTolerance',1e-9,'IntegerTolerance',1e-6);

Храните историю вычисленных истинных квадратичных и слабых переменных для печати.

history = [truequadratic,zslack];

Вычислите квадратичные и слабые значения. Если они отличаются, добавьте еще одно линейное ограничение и решите снова.

В синтаксисе панели инструментов каждая новая линейная зависимость Ax≤b происходит от линейной аппроксимации

- xkTQxk+2xkTQx-z≤0.

Вы видите, что новая строка A = 2xkTQ и новый элемент в b = xkTQxk, где член z представлен коэффициентом -1 в A.

После поиска нового решения используйте линейное ограничение на полпути между старым и новым решениями. Этот эвристический способ включения линейных ограничений может быть более быстрым, чем простой выбор нового решения. Чтобы использовать решение вместо эвристики половины пути, прокомментируйте строку «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)')

Figure contains an axes. The axes with title Quadratic and linear approximation (slack) contains 2 objects of type line. These objects represent Quadratic, Slack.

Каково качество решения MILP? output структура содержит эту информацию. Проверьте абсолютный разрыв между внутренними вычисленными границами цели в решении.

disp(output.absolutegap)
     0

Абсолютный промежуток равен нулю, что указывает на точность решения MILP.

Постройте график оптимального распределения. Использовать xLinInt(xvars), не assets, потому что assets может не удовлетворять ограничениям при использовании промежуточного обновления.

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

Figure contains an axes. The axes with title Optimal Asset Allocation contains an object of type bar.

Легко увидеть, что все ненулевые соотнесения основных средств находятся между полунепрерывными границами 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 для непосредственной обработки полунепрерывных ограничений и ограничений количества элементов см. в разделе Оптимизация портфеля с полунепрерывными ограничениями и ограничениями количества элементов (Financial Toolbox).

Связанные темы