www.gusucode.com > optim 案例源码 matlab代码程序 > optim/PortfolioMIQPExample.m

    %% Mixed-Integer Quadratic Programming Portfolio Optimization
%
% This example shows how to solve a Mixed-Integer Quadratic Programming
% (MIQP) portfolio optimization problem using the |intlinprog|
% Mixed-Integer Linear Programming (MILP) solver. The idea is to
% iteratively solve a sequence of MILP problems that locally approximate
% the MIQP problem.

% Copyright 2015 The MathWorks, Inc.
%% Problem Outline
%
% As Markowitz showed ("Portfolio Selection," J. Finance Volume 7, Issue 1,
% pp. 77-91, March 1952), you can express many portfolio optimization
% problems as quadratic programming problems. Suppose that you have a set
% of |N| assets and want to choose a portfolio, with $x(i)$ being the
% fraction of your investment that is in asset $i$. If you know the vector
% $r$ of mean returns of each asset, and the covariance matrix $Q$ of the
% returns, then for a given level of risk-aversion $\lambda$ you maximize
% the risk-adjusted expected return:
%
% $$\max_x(r^Tx - \lambda x^TQx).$$
%
% The |quadprog| solver addresses this quadratic programming problem.
% However, in addition to the plain quadratic programming problem, you
% might want to restrict a portfolio in a variety of ways, such as:
%
% * Having no more than |M| assets in the portfolio, where |M <= N|.
% * Having at least |m| assets in the portfolio, where |0 < m <= M|.
% * Having _semicontinuous_ constraints, meaning either $x(i) = 0$, or
% $f{\mathrm{min}}\le x(i)\le f{\mathrm{max}}$ for some fixed fractions
% $f{\mathrm{min}} > 0$ and $f{\mathrm{max}} \ge f{\mathrm{min}}$.
%
% You cannot include these constraints in |quadprog|. The difficulty is the
% discrete nature of the constraints. Furthermore, while the mixed-integer
% linear programming solver |intlinprog| does handle discrete constraints,
% it does not address quadratic objective functions.
%
% This example constructs a sequence of MILP problems that satisfy the
% constraints, and that increasingly approximate the quadratic objective
% function. While this technique works for this example, it might not apply
% to different problem or constraint types.
%
% Begin by modeling the constraints.
%% Modeling Discrete Constraints
%
% $x$ is the vector of asset allocation fractions, with $0\le x(i)\le 1$
% for each $i$. To model the number of assets in the portfolio, you need
% 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$ vector to be a binary variable, and impose the linear constraints
%
% $$v(i)f{\mathrm{min}} \le x(i)\le v(i)f{\mathrm{max}}.$$
%
% These inequalities both enforce that $x(i)$ and $v(i)$ are zero at
% exactly the same time, and they also enforce that $f{\mathrm{min}} \le
% x(i)\le f{\mathrm{max}}$ whenever $x(i) > 0$.
%
% Also, to enforce the constraints on the number of assets in the
% portfolio, impose the linear constraints
%
% $$m\le \sum_i v(i)\le M.$$
%% Objective and Successive Linear Approximations
%
% As first formulated, you try to maximize the objective function. However,
% all Optimization Toolbox(TM) solvers minimize. So formulate the problem
% as minimizing the negative of the objective:
%
% $$ \min_x \lambda x^TQx - r^Tx .$$
%
% This objective function is nonlinear. 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\mbox{ such that } x^TQx - z \le 0,\ z\ge 0 .$$
%
% As you iteratively solve MILP approximations, you include new linear
% constraints, each of which approximates the nonlinear constraint locally
% near the current point. In particular, for $x = x_0 + \delta$ where $x_0$
% is a constant vector and $\delta$ is a variable vector, the first-order
% Taylor approximation to the constraint is
%
% $$x^TQx - z = x_0^TQx_0 + 2x_0^TQ\delta - z + O(|\delta|^2).$$
%
% Replacing $\delta$ by $x - x_0$ gives
%
% $$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 \le 0.$$
%
% This has the form $Ax\le b$, where $A = 2x_k^TQ$, there is a $-1$
% multiplier for the $z$ term, and $b = x_k^TQx_k$.
%
% This method of adding new linear constraints to the problem is called a
% cutting plane method. For details, see 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) Problem Formulation
%
% To express problems for the |intlinprog| solver, you need to do the
% following:
%
% * Decide what your variables represent
% * Express lower and upper bounds in terms of these variables
% * Give linear equality and inequality matrices
%
% Have the first $N$ variables represent the $x$ vector, the next $N$
% variables represent the binary $v$ vector, and the final variable
% represent the $z$ slack variable. There are $2N + 1$ variables in the
% problem.
%
% Load the data for the problem. This data has 225 expected returns in the
% vector |r| and the covariance of the returns in the 225-by-225 matrix
% |Q|. The data is the same as in the Using Quadratic Programming on
% Portfolio Optimization Problems example.
load port5
r = mean_return;
Q = Correlation .* (stdDev_return * stdDev_return');
%%
% Set the number of assets as |N|.
N = length(r);
%%
% Set indexes for the variables
xvars = 1:N;
vvars = N+1:2*N;
zvar = 2*N+1;
%%
% The lower bounds of all the |2N+1| variables in the problem are zero. The
% upper bounds of the first |2N| variables are one, and the last variable
% has no upper bound.
lb = zeros(2*N+1,1);
ub = ones(2*N+1,1);
ub(zvar) = Inf;
%%
% Set the number of assets in the solution to be between 100 and 150.
% Incorporate this constraint into the problem in the form, namely
%
% $$m\le \sum_i v(i)\le M,$$
%
% by writing two linear constraints of the form $Ax \le b$:
%
% $$\sum_i v(i) \le M$$
%
% $$\sum_i -v(i) \le -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;
%%
% Include semicontinuous constraints. Take the minimal nonzero fraction of
% assets to be |0.001| for each asset type, and the maximal fraction to be
% |0.05|.
fmin = 0.001;
fmax = 0.05;
%%
% Include the inequalities $x(i)\le f{\mathrm{max}}(i)*v(i)$ and
% $f{\mathrm{min}}(i)*v(i)\le x(i)$ as linear inequalities.
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)];
%%
% Include the constraint that the portfolio is 100% invested, meaning $\sum
% x_i = 1$.
Aeq = zeros(1,2*N+1); % Allocate Aeq matrix
Aeq(xvars) = 1;
beq = 1;
%%
% Set the risk-aversion coefficient $\lambda$ to |100|.
lambda = 100;
%%
% Define the objective function $\lambda z - r^Tx$ as a vector. Include
% zeros for the multipliers of the $v$ variables.
f = [-r;zeros(N,1);lambda];
%% Solve the Problem
% To solve the problem iteratively, begin by solving the problem with the
% current constraints, which do not yet reflect any linearization. The
% integer constraints are in the |vvars| vector.
options = optimoptions(@intlinprog,'Display','off'); % Suppress iterative display
[xLinInt,fval,exitFlagInt,output] = intlinprog(f,vvars,A,b,Aeq,beq,lb,ub,options);
%%
% Prepare a stopping condition for the iterations: stop when the slack
% variable $z$ is within 0.01% of the true quadratic value.
thediff = 1e-4;
iter = 1; % iteration counter
assets = xLinInt(xvars); % the x variables
truequadratic = assets'*Q*assets;
zslack = xLinInt(zvar); % slack variable value
%%
% Keep a history of the computed true quadratic and slack variables for
% plotting.
history = [truequadratic,zslack];
%%
% Compute the quadratic and slack values. If they differ, then add another
% linear constraint and solve again.
%
% In toolbox syntax, each new linear constraint $Ax\le b$ comes from the
% linear approximation
%
% $$-x_k^TQx_k + 2x_k^TQx - z \le 0.$$
%
% You see that the new row of $A = 2x_k^TQ$ and the new element in  $b =
% x_k^TQx_k$, with the $z$ term represented by a -1 coefficient in $A$.
%
% After you find a new solution, use a linear constraint halfway between
% the old and new solutions. This heuristic way of including linear
% constraints can be faster than simply taking the new solution. To use the
% solution instead of the halfway heuristic, comment the "Midway" line
% below, and uncomment the following one.
while abs((zslack - truequadratic)/truequadratic) > thediff % relative error
    newArow = horzcat(2*assets'*Q,zeros(1,N),-1); % Linearized constraint
    A = [A;newArow];
    b = [b;truequadratic];
    % 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 = assets'*Q*assets;
    zslack = xLinInt(zvar);
    history = [history;truequadratic,zslack];
    iter = iter + 1;
end
%% Examine the Solution and Convergence Rate
% Plot the history of the slack variable and the quadratic part of the
% objective function to see how they converged.
plot(history)
legend('Quadratic','Slack')
xlabel('Iteration number')
title('Quadratic and linear approximation (slack)')
%%
% What is the quality of the MILP solution? The |output| structure contains
% that information. Examine the absolute gap between the
% internally-calculated bounds on the objective at the solution.
disp(output.absolutegap)
%%
% The absolute gap is zero, indicating that the MILP solution is accurate.
%%
% Plot the optimal allocation. Use |xLinInt(xvars)|, not |assets|, because
% |assets| might not satisfy the constraints when using the midway update.
bar(xLinInt(xvars))
grid on
xlabel('Asset index')
ylabel('Proportion of investment')
title('Optimal asset allocation')
%%
% You can easily see that all nonzero asset allocations are between the
% semicontinuous bounds $f{\mathrm{min}} = 0.001$ and $f{\mathrm{max}} =
% 0.05$.
%
% How many nonzero assets are there? The constraint is that there are
% between 100 and 150 nonzero assets.
sum(xLinInt(vvars))
%%
% What is the expected return for this allocation, and the value of the
% risk-adjusted return?
fprintf('The expected return is %g, and the risk-adjusted return is %g.\n',...
    r'*xLinInt(xvars),-fval)

%%
% More elaborate analyses are possible by using features specifically
% designed for portfolio optimization in Financial Toolbox(TM).