## Optimization Toolbox |

This example shows how to solve portfolio optimization problems using the interior-point quadratic programming algorithm in `quadprog`

. The function `quadprog`

belongs to Optimization Toolbox™.

The matrices that define the problems in this example are dense; however, the interior-point algorithm in `quadprog`

can also exploit sparsity in the problem matrices for increased speed. See the examples section of Optimization Toolbox documentation for a sparse example.

Suppose that there are different assets. The rate of return of asset is a random variable with expected value . The problem is to find what fraction to invest in each asset in order to minimize risk, subject to a specified minimum expected rate of return.

Let denote the covariance matrix of rates of asset returns.

The classical mean-variance model consists of minimizing portfolio risk, as measured by

subject to a set of constraints.

The expected return should be no less than a minimal rate of portfolio return that the investor desires,

the sum of the investment fractions 's should add up to a total of one,

and, being fractions (or percentages), they should be numbers between zero and one,

Since the objective to minimize portfolio risk is quadratic, and the constraints are linear, the resulting optimization problem is a quadratic program, or QP.

Let us now solve the QP with 225 assets. The dataset is from the OR-Library [Chang, T.-J., Meade, N., Beasley, J.E. and Sharaiha, Y.M., "Heuristics for cardinality constrained portfolio optimisation" Computers & Operations Research 27 (2000) 1271-1302].

We load the dataset and then set up the constraints in a format expected by `quadprog`

. In this dataset the rates of return
range between -0.008489 and 0.003971; we pick a desired return
in between, e.g., 0.002 (0.2 percent).

% Load dataset stored in a MAT-file. load('port5.mat','Correlation','stdDev_return','mean_return') % Calculate covariance matrix from correlation matrix. Covariance = Correlation .* (stdDev_return * stdDev_return'); nAssets = numel(mean_return); r = 0.002; % number of assets and desired return Aeq = ones(1,nAssets); beq = 1; % equality Aeq*x = beq Aineq = -mean_return'; bineq = -r; % inequality Aineq*x <= bineq lb = zeros(nAssets,1); ub = ones(nAssets,1); % bounds lb <= x <= ub c = zeros(nAssets,1); % objective has no linear term; set it to zero

In order solve the QP using the interior-point algorithm, we set the option Algorithm to 'interior-point-convex'.

options = optimoptions('quadprog','Algorithm','interior-point-convex');

We now set some additional options, and call the solver quadprog.

% Set additional options: turn on iterative display, and set a tighter optimality termination tolerance. options = optimoptions(options,'Display','iter','TolFun',1e-10); % Call solver and measure wall-clock time. tic [x1,fval1] = quadprog(Covariance,c,Aineq,bineq,Aeq,beq,lb,ub,[],options); toc % Plot results. plotPortfDemoStandardModel(x1)

First-order Total relative Iter f(x) Feasibility optimality error 0 7.212813e+00 1.227e+02 1.196e+00 6.769e+02 1 8.160874e-04 3.615e-01 3.522e-03 4.321e+00 2 7.220766e-04 3.593e-01 3.500e-03 3.247e+00 3 4.309434e-04 9.991e-02 9.734e-04 1.660e+00 4 4.734300e-04 2.220e-16 4.242e-06 3.823e-01 5 4.719034e-04 4.441e-16 8.003e-07 5.754e-03 6 3.587475e-04 4.441e-16 3.677e-07 3.650e-03 7 3.131814e-04 2.220e-16 9.587e-08 1.197e-03 8 2.760174e-04 4.441e-16 1.521e-08 2.652e-04 9 2.345751e-04 2.220e-16 4.110e-09 1.161e-04 10 2.042487e-04 2.220e-16 6.423e-09 2.748e-05 11 1.961775e-04 2.220e-16 6.068e-10 3.121e-06 12 1.949281e-04 4.441e-16 4.280e-12 2.705e-08 Minimum found that satisfies the constraints. Optimization completed because the objective function is non-decreasing in feasible directions, to within the selected value of the function tolerance, and constraints are satisfied to within the default value of the constraint tolerance. Elapsed time is 0.259245 seconds.

We now add to the model group constraints that require that 30% of the investor's money has to be invested in assets 1 to 75, 30% in assets 76 to 150, and 30% in assets 151 to 225. Each of these groups of assets could be, for instance, different industries such as technology, automotive, and pharmaceutical. The constraints that capture this new requirement are

% Add group constraints to existing equalities. Groups = blkdiag(ones(1,nAssets/3),ones(1,nAssets/3),ones(1,nAssets/3)); Aineq = [Aineq; -Groups]; % convert to <= constraint bineq = [bineq; -0.3*ones(3,1)]; % by changing signs % Call solver and measure wall-clock time. tic [x2,fval2] = quadprog(Covariance,c,Aineq,bineq,Aeq,beq,lb,ub,[],options); toc % Plot results, superimposed to results from previous problem. plotPortfDemoGroupModel(x1,x2);

First-order Total relative Iter f(x) Feasibility optimality error 0 7.212813e+00 1.227e+02 3.540e-01 6.737e+02 1 7.004556e-03 2.901e+00 8.367e-03 1.814e+01 2 9.181962e-04 4.096e-01 1.181e-03 4.220e+00 3 7.515047e-04 3.568e-01 1.029e-03 5.242e+00 4 4.238346e-04 9.006e-02 2.597e-04 1.510e+00 5 3.695008e-04 1.910e-04 1.342e-05 3.024e-01 6 3.691407e-04 6.146e-07 6.817e-08 9.741e-04 7 3.010636e-04 7.692e-08 1.837e-08 4.322e-04 8 2.669065e-04 1.088e-08 5.475e-09 2.007e-04 9 2.195767e-04 8.123e-10 2.814e-08 8.957e-05 10 2.102910e-04 2.840e-10 1.037e-08 3.391e-05 11 2.060985e-04 6.714e-11 2.877e-09 1.518e-05 12 2.015107e-04 2.220e-16 1.522e-10 1.389e-06 13 2.009670e-04 6.661e-16 5.264e-13 4.771e-09 Minimum found that satisfies the constraints. Optimization completed because the objective function is non-decreasing in feasible directions, to within the selected value of the function tolerance, and constraints are satisfied to within the default value of the constraint tolerance. Elapsed time is 0.248330 seconds.

We see from the second bar plot that, as a result of the additional group constraints, the portfolio is now more evenly distributed across the three asset groups than the first portfolio. This imposed diversification also resulted in a slight increase in the risk, as measured by the objective function (see column labelled "f(x)" for the last iteration in the iterative display for both runs).

In order to show how `quadprog`

's interior-point algorithm behaves on a larger problem, we'll use a 1000-asset randomly generated dataset. We generate a random correlation matrix (symmetric, positive-semidefinite, with ones on the diagonal) using the `gallery`

function in MATLAB®.

% Reset random stream for reproducibility. rng(0,'twister'); nAssets = 1000; % desired number of assets % Generate means of returns between -0.1 and 0.4. a = -0.1; b = 0.4; mean_return = a + (b-a).*rand(nAssets,1); % Generate standard deviations of returns between 0.08 and 0.6. a = 0.08; b = 0.6; stdDev_return = a + (b-a).*rand(nAssets,1); % Correlation matrix, generated using Correlation = gallery('randcorr',nAssets). % (Generating a correlation matrix of this size takes a while, so we load % a pre-generated one instead.) load('correlationMatrixDemo.mat','Correlation'); % Calculate covariance matrix from correlation matrix. Covariance = Correlation .* (stdDev_return * stdDev_return');

We now define the standard QP problem (no group constraints here) and solve.

r = 0.15; % desired return Aeq = ones(1,nAssets); beq = 1; % equality Aeq*x = beq Aineq = -mean_return'; bineq = -r; % inequality Aineq*x <= bineq lb = zeros(nAssets,1); ub = ones(nAssets,1); % bounds lb <= x <= ub c = zeros(nAssets,1); % objective has no linear term; set it to zero % Call solver and measure wall-clock time. tic x3 = quadprog(Covariance,c,Aineq,bineq,Aeq,beq,lb,ub,[],options); toc

First-order Total relative Iter f(x) Feasibility optimality error 0 2.142849e+01 5.490e+02 3.032e+00 1.102e+04 1 9.378552e-03 6.439e+00 3.556e-02 3.143e+02 2 1.128129e-04 3.706e-03 2.047e-05 9.457e+00 3 1.118804e-04 1.853e-06 1.171e-07 5.260e-03 4 8.490176e-05 7.650e-08 7.049e-09 3.387e-04 5 3.364597e-05 2.220e-16 1.037e-09 1.641e-04 6 1.980189e-05 2.220e-16 8.466e-11 3.748e-05 Minimum found that satisfies the constraints. Optimization completed because the objective function is non-decreasing in feasible directions, to within the selected value of the function tolerance, and constraints are satisfied to within the default value of the constraint tolerance. Elapsed time is 1.813593 seconds.

This example illustrates how to use the interior-point algorithm in `quadprog`

on a portfolio optimization problem, and shows the algorithm running times on quadratic problems of different sizes. The runs were made on a 64-bit, 8-core, 3 GHz Intel® Xeon® cpu with 12 gigabytes of memory running Linux® operating system.

More elaborate analyses are possible by using features specifically designed for portfolio optimization in Financial Toolbox™.