【MILP】Mixed-Integer Quadratic Programming portfolio optimzation

本文链接：https://blog.csdn.net/qq_18822147/article/details/114586747

本文介绍了一个经典的Markowitz资产分配问题，并提出了一种通过引入松弛变量和使用整数线性规划求解的方法。通过线性逼近非线性约束，迭代地解决混合整数线性规划(MILP)问题。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

Navigator

Problem
Objective and Successive linear Approximations
Code
cvx code
Reference

Problem

The classical Markowitz problem is as follows:
$\max_x r^Tx-\lambda x^TQx$
The quadprog solver addresses this quadratic programming problem. Some cardinality constraints would be added. $x$ is the vector of asset allocation fractions, with $0\leq x(i)\leq 1$ for each $i$ . The indicator variables $v$ such that $v (i) = 0$ when $x (i) = 0$ , and $v (i) = 1$ when $x (i) > 0$ . To get variables that satisfy this restriction, set the $v$ vectror to be a binary variable, and impose the linear constraints.
$v(i)f_{\min}\leq x(i)\leq v(i)f_{\max}$
and
$m\leq \sum v(i)\leq M$

Objective and Successive linear Approximations

The objective function is
$\min_x\lambda x^TQx-r^Tx$
However, the intlinprog MILP solver requires a linear objective function. There is a standard technique to reformulate this problem into one with linear objective and nonlinear constraints. Introduce a slack variable $z$ to represent the quadratic term.
$\min_{x, z} \lambda z-r^Tx\\ \begin{cases} x^TQx-z\leq 0\\ z\geq 0 \end{cases}$
To solve MILP approximations iteratively, we could approximate the nonlinear constraint locally near the current point
$x^TQx-z=x_0^TQx_0+2x_0^TQ\delta-z+O(|\delta|^2)$
Replacing $\delta$ by $x-x_0$
$x^TQx-z = -x_0^TQx_0+2x_0^TQx-z+O(|x-x_0|^2)$
For each intermediate solution $x_k$ you introduce a new linear constraint in $x$ and $z$ as the linear part of the expression above:
$-x_k^TQx_k+2x_k^TQx-z\leq 0$

Code

%% integer programming
clc;
clear all;
close all;

%%
load port5;
r = mean_return;
Q = Correlation .* (stdDev_return * stdDev_return');
N = length(r);
% indexes
xvars = 1:N;
vvars = N+1:2*N;
zvar = 2*N+1;

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

% set the number of assets in the solution to between 100 and 150.
M = 150;
m = 100;
A = zeros(1, 2*N+1);
A(vvars) = 1; % A*x represents the sum of v(i)
A = [A; -A];
b = zeros(2, 1);
b(1) = M;
b(2) = -m;

% set the minimal nonzero fraction of assets to be 0.05 and the maximal to
% be 0.20
fmin = 0.001;
fmax = 0.05;

% inequalities fmin(i)*v(i)<=x(i)<= fmax(i)*v(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)];

% constraint sum(x_i)=1
Aeq = zeros(1, 2*N+1);
Aeq(xvars) = 1;
beq = 1;

% risk-aversion coefficient lambda
lambda = 100;

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

%% solve
ops = optimoptions(@intlinprog, 'Display', 'off');
[xLinInt, fval, exitFlag, output] = intlinprog(f, vvars, A, b, Aeq, beq, lb, ub, ops);

%% stopping condition
thediff = 1e-4;
iter = 1;
assets = xLinInt(xvars); 
trueVar = assets'*Q*assets;
zslack = xLinInt(zvar);

cvx code

%% MLIP
n = 10;
k_high = 8;
k_low = 2;
fmax = 0.30;
fmin = 0.05;
mu = abs(randn(n, 1));
Sigma = randn(n, n)/10;
Sigma = Sigma'*Sigma;

cvx_solver mosek
cvx_begin
    variable w(n);
    variable bi(n) binary;
    ret = mu'*w;
    risk = quad_form(w, Sigma);
    maximize(ret-risk);
    subject to 
        ones(1, n)*w==1;
        w>=0;
        w<=bi;
        % k_high>=nnz(w)>=k_low;
        % norm(w, Inf)<=fmax;
        fmin*bi<=w;
        fmax*bi>=w;
        sum(bi)>=k_low;
        sum(bi)<=k_high;
cvx_end

w=value(w);