Find matrix that maximizes formula output.

rates = [ 0.0174 0.0188 0.0186 0.0173 0.0192 0.0179 0.0192 0.0185 0.0175 0.0196 0.0184 0.0195 0.0187 0.0178 0.0199 0.0187 0.0199 0.0189 0.0182 0.0203 0.0188 0.0204 0.0191 0.0186 0.0208 0.0190 0.0210 0.0194 0.0189 0.0212 0.0191 0.0215 0.0196 0.0193 0.0216 0.0192 0.0217 0.0196 0.0196 0.0219 0.0193 0.0220 0.0197 0.0200 0.0220 0.0194 0.0000 0.0198 0.0203 0.0223 0.0194 0.0000 0.0199 0.0204 0.0225 0.0195 0.0000 0.0200 0.0207 0.0226];

weights = [ 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000];

yield = sum(sum(rates.* weights)) * 100;

It needs to change the weights to an optimal array that maximizes the yield formula. I've tried the fmin functions, I've tried the linprog function, I attempted simulink... just don't know how to make it use a matrix as the changing variable.

EX: in Excel, solver will return this matrix to maximize yield.

weights = [ 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.2675 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.3000];

Answer 1

[Walter Roberson]

Assuming this is something like an asset mix where you want the total weight to be 1, then:

weights = linprog(-rates.', [], [], ones(1,60), 1, zeros(1,60), ones(1,60))

This returns a vector that is all 0 except for the last entry, corresponding to rate 0.0226, which is the maximum entry. This reflects the fact that you would get the highest return by concentrating everything on the stock with the highest rate.

If you put a different upper bound, limiting the exposure to any one asset, such as ones(1,60)/20 then the result you get out will be equivalent to taking the maximum exposure in the 20 largest-yielding assets.

[Damani Brown]

You are a genius. Yes, the Excel model has weights on asset classes. The matrix consists of different assets and was going to be what I incorporated next.

Thanks for this, I can now finish my optimizer.

[Walter Roberson]

If you have minimum and maximum exposure to assets, then you would code that information in place of zeros(1,60), ones(1,60) .

But notice that as long as the purchases do not affect each other, then the maximum is always going to be found by starting with the minimum exposure permitted on all of the items, and then going through and raising the exposure to maximum on all of the remaining items from the best yielding down towards the less good, until eventually the sum totals 1. No real optimization is required, just initialization, a sort, and a for loop.

[Damani Brown]

Yes, we have the solver working as intended in Excel with constraints. However, before completion the number of changing variables were well over 400, greater than Excel's 200 max. I don't believe Matlab has such a constraint and we will utilize it's full power.

Currently this solver will optimize exposure to bond sectors: [Agency, Treasury Bill, Treasury Note, CD/YCD, Corporate]

with constraints: [30%, 100%, 100%, 40%, 20%]

Given the rate matrix, where should it put the money with given the constraints.

For a full view of what I am doing... the matrix is set up as such.

Sector[agency, tbill, tnote, cds, corp]

[6 days[rate, rate, rate, rate, rate]

[30 days[rate, rate, rate, rate, rate]

So every row are rates given a different maturity bucket. every column is a different security class.

the sum of each column cannot exceed the constraint. Ex: Agency bonds cannot exceed 30% of the fund balance. Assume we are investing 100% of the fund today. It can invest up to the constraints in any sector -- it can even invest the entire fund in tbills or tnotes if it chooses so since exposure can be 100%.

That's just an overview of what I am doing here just incase the rate table is a little confusing.

[Damani Brown]

I need to add a constraint to where the sum of the weight matrix does not exceed constraints = [0.4000 0.1000 1.0000 1.0000 0.3000];

each column(s) sum cannot exceed the 1x5 matrix. I'll research to see if this is possible.

[Walter Roberson]

For "cannot exceed", use the A and b entries. For example,

A = [zeros(1,80*0), ones(1,80), zeros(1,400-80*1);
     zeros(1,80*1), ones(1,80), zeros(1,400-80*2);
     zeros(1,80*2), ones(1,80), zeros(1,400-80*3);
     zeros(1,80*3), ones(1,80), zeros(1,400-80*4);
     zeros(1,80*4), ones(1,80), zeros(1,400-80*5)];
b = [0.4000 0.1000 1.0000 1.0000 0.3000];

[Damani Brown]

The number of columns in A must be the same as the number of elements of f which is our rates matrix converted into a 1 x 60 matrix.

[Walter Roberson]

So change the 80's to 12's.

[Damani Brown]

You are a genius. I also needed to convert the 400 to 60. Thanks, worked great.

[Walter Roberson]

Use variables to represent the slice width and the total size: it makes the code easier to adapt.

As you are using linprog() you might be interested in the newer Problem Based approach; https://www.mathworks.com/help/optim/problem-based-lp-milp.html

[Damani Brown]

You're right, I stumbled across a guy who was attempting something similar to me. https://www.mathworks.com/matlabcentral/answers/365929-how-do-i-solve-for-portfolio-weights-to-achieve-specific-min-max-objective?s_tid=srchtitle

[Damani Brown]

Walter, since you're a genius. Is there a name for the matrix that looks like an non-square identity matrix. I've seen this before in Excel.

A = [1 1 0 0 0 0; 0 0 1 1 0 0; 0 0 0 0 1 1]

[Walter Roberson]

I do not know of any particular name for that, but

kron(eye(3),ones(1,2))

can reproduce it.

[Damani Brown]

That's a nifty trick. Thanks!

[Damani Brown]

We have named the matrix: The Walter Matrix.

Answer 2

[Damani Brown]

Tried to convert it to the new optimization method, got close until I attempted to add the second constraint.

R = [
0.0174  0.0188  0.0186  0.0172  0.0192
0.0179  0.0192  0.0185  0.0175  0.0196
0.0184  0.0195  0.0187  0.0178  0.0199
0.0187  0.0199  0.0189  0.0182  0.0203
0.0188  0.0204  0.0191  0.0186  0.0208
0.0190  0.0210  0.0194  0.0189  0.0212
0.0191  0.0215  0.0196  0.0193  0.0216
0.0192  0.0217  0.0196  0.0196  0.0219
0.0193  0.0220  0.0197  0.0200  0.0220
0.0194  0.0000  0.0198  0.0203  0.0223
0.0194  0.0000  0.0199  0.0204  0.0225
0.0195  0.0000  0.0200  0.0207  0.0226];

b = [.1000 .4000 1 1 .3000];
R = reshape(R,60,[]);
walter_matrix = kron(eye(5),ones(1,12));
weights = optimvar('weights', 60, 1, 'UpperBound', 1, 'LowerBound', 0);

prob = optimproblem('ObjectiveSense','max');
prob.Objective = sum(R .* weights);
prob.Constraints.Constr1= sum(weights) == 1;
prob.Constraints.Constr2 = walter_matrix * weights <= b;
sol = solve(prob);

Returned: Error using <= Argument dimensions 5-by-1 and 1-by-5 must agree.

Error in solve (line 23) prob.Constraints.Constr2 = walter_matrix * weights <= b;

If I used the ".*" method, it will return Error using optim.internal.problemdef.ElementwiseOperator/checkIsValid Argument dimensions 5-by-60 and 60-by-1 must agree.

[Damani Brown]

I figured out the answer.

[Damani Brown]

Posting OptimVar solution for anyone trying to maximize return on a portfolio.

% Days to maturity
days = [
    6.00
    37.00
    67.00
    98.00
    128.00
    159.00
    190.00
    218.00
    249.00
    279.00
    310.00
    340.00];

% Rates per sector
rates = [
% AGD   CP      T-BILL  T-NOTE  YCD  
0.0174  0.0188  0.0186  0.0172  0.0192
0.0179  0.0192  0.0185  0.0175  0.0196
0.0184  0.0195  0.0187  0.0178  0.0199
0.0187  0.0199  0.0189  0.0182  0.0203
0.0188  0.0204  0.0191  0.0186  0.0208
0.0190  0.0210  0.0194  0.0189  0.0212
0.0191  0.0215  0.0196  0.0193  0.0216
0.0192  0.0217  0.0196  0.0196  0.0219
0.0193  0.0220  0.0197  0.0200  0.0220
0.0194  0.0000  0.0198  0.0203  0.0223
0.0194  0.0000  0.0199  0.0204  0.0225
0.0195  0.0000  0.0200  0.0207  0.0226];

% Walter Matrix
walter_matrix = kron(eye(5),ones(1,12));

% Constraints
agd = 0.1000;
cp = 0.4000;
ustb = 1.0000;
ustn = 1.0000;
ycd = 0.3000;
constraints = [agd cp ustb ustn ycd];

% Optimize Portfolio
rates = reshape(rates,60,[]); % convert to 60 x 1 matrix
weights = optimvar('weights',60,1,'LowerBound',0,'UpperBound',1);
port = optimproblem('ObjectiveSense','max');
port.Objective = sum(rates .* weights);
port.Constraints.total = sum(weights) == 1;
port.Constraints.ips = walter_matrix * weights <= constraints.';
solution = solve(port);

% Print solution weight matrix
disp([days reshape(solution.weights,12,[])]);

% Compute Yield
yield = sum(rates .* solution.weights) * 100;
fprintf("Optimization yield is %0.2f%%\n",yield);