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).