Linear programming and cash flow matching

I am trying to solve a simple large scale linear /non-linear programming problem to match liability cash flows with asset cash flows. I am thinking between MATLAB and R to solve this problem. I have explained the problem below and I was wondering whether MATLAB can be used to solve this.

Can anyone please let me know if and how MATLAB can be used to solve this problem? If any particular MATLAB optimization toolboxes will be needed then please let me know.

A fixed set of liability cash flows are given, say for the next 20 years. We need to select appropriate bonds to match these liability cash flows, ignoring interest rate risk and other market risks. We are given market value of bonds and their cash flows based on which we need to decide which bonds to select. The objective function is to minimise value of liabilities, which is calculated by discounting liability cash flows at a flat discount rate, say 1% + a spread. This spread depends upon internal rate of return of the assets portfolio, e.g. spread = 50% of internal rate of return on the asset portfolio (the higher the spread, the lower the value of liabilities. Constraints are: weights of bonds can be either zero or 1 (i.e. either we invest in a bond or we don’t) and asset cash flows can be lower than liability cash flows in certain years but cumulatively say at every 3 years asset cash flows should exceed liability cash flows.

Answer 1

[Alan Weiss]

I don't know for sure, but perhaps this example is relevant.

Alan Weiss

MATLAB mathematical toolbox documentation

[L Smith]

Thank you Alan for your help. I have found following example from a webinar which is quite relevant for the type of problem I am trying to solve: https://uk.mathworks.com/matlabcentral/fileexchange/65158-linear-and-mixed-integer-linear-programming-in-matlab?s_cid=ME_prod_FX

In the example, there were 5 different bonds to choose from. Can you please let me know whether this approach can be use in a large scale problem? For example, if there are 1000 different types of bonds to choose from?

Generally, are MATLAB optimization tools powerful enough for large scale problems? I have heard mixed reviews from some people with some commenting that MATLAB is not suitable for large-scale problem but it seems like there are a number of optimization toolboxes available from Mathworks.

[Torsten]

For large-scale problems, usually CPLEX is the solver of choice.

[Mary Fenelon]

Are you interested in original formulation (linear programming, LP) or the round lots formulation (mixed-integer programming, MILP)?

It's easier to make predictions about LP solution times; MATLAB will solve a problem with 1000 bonds quickly. Predictions about MILP solution times are much harder. For a given formulation, it might be easy to solve with one set of data and difficult with another.

That said, this particular problem is easy to solve, either the LP or MILP formulation. I generated random data for prices and durations for 8 periods and 1000 bonds using this code:

nt = 8;
nb = 1000;

prices = round(100*rand(nb,1),2);
durations = randi(nt-2,[nb,1])+2;
coupons = round((2 + 2*rand(nb,1)).*prices/100,2)';

cashFlows = repmat(coupons,nt,1);
for i = 1:nb
    cashFlows(durations(i),i) = ceil(prices(i)+coupons(i));
    cashFlows(durations(i)+1:nt,i) = 0;
end

For a few different random seeds, on my laptop the LP always solved in less than 0.05 seconds; the MILP solve time varied between about 0.2 and 40 seconds.

You may have some additional constraints on which bonds to use; these may not have much effect on the LP formulation but could make the MILP formulation easier or harder. Do a web search on "integer programming peformance variability" to find some references on why it's so hard to predict.

[L Smith]

Thanks Mary. I am interested in the mixed-integer programming example. The problem that I am trying to solver is slightly different from the example covered in the webinar. Here is a summary of the problem

• Objective function is to minimize the value of liabilities- calculated as value of obligation cash flows discounted at the internal rate of return on the bonds portfolio selected.
• Optimization variables (i.e. bonds) are binary-we can only have one of each as max
• Constraints are that total cash flows of bonds selected should be greater than or equal to obigation cash flows

However, I am running into a problem with calculating internal rate of return on bonds within the optimization process. The code below is adapted from your code in the webinar and the line at which MATLAB gives an error is in bold font.

Q1. Can you please let me know if there is a way to fix this? The problem seems to be with the irr function which doesn't accept any arguments other than static cash flows.

Q2. Given this objective function and the constraints explained above, do you think MATLAB optimization will be fast for a large scale problem, say if number of bonds are 1000?

clear
% Define data 
prices = [99.74 91.22 98.71 103.75 97.15];
cashFlows = [4 5 2.5 5 4; 4 5 2.5 5 4;4 5 2.5 5 4;4 5 2.5 5 4;4 5 102.5 5 4;4 5 0 105 104;
4 105 0 0 0;104 0 0 0 0];
obligations = [5 7 7 6 8 7 20 0]';
nt=size(cashFlows,1)
nb=size(cashFlows,2)
matchprob=optimproblem;
bonds = optimvar('bonds',nb,'Type','integer','LowerBound',0,'UpperBound',1);
matchprob.ObjectiveSense = 'minimize';
%Internal rate of return calculation 
cashFlows2 = [-prices; cashFlows]
f = [cashFlows2*bonds]'
LiabDiscountR=irr(f);   
matchprob.Objective = pvar(obligations,LiabDiscountR);
matchprob.Constraints.cover = cashFlows*bonds >= obligations;
showproblem(matchprob)
milpSol = solve(matchprob);
milpSol.bonds

[Walter Roberson]

irr is defined on a vector or matrix, not on on OptimizationExpression

[L Smith]

Do you mean I cannot use irr on Optimization expression? The Optimization objective is essentially to maximise internal rate of return. Can you please let me know how I can correct the code above?

[Walter Roberson]

I do mean that you cannot use irr() on an OptimizationExpression variable. irr() was designed long long before OptimizationExpressions existed.

As far as I know, there is no way to apply an arbitrary function to an OptimizationExpression

If it happened to be the case that irr() was what you were trying to optimize subject to a bunch of constraints, then you could potentially construct everything else using optimization expressions, and then have MATLAB convert the optimzation setup into an optimization Problem -- and then wrap an irr() call around the optimization function call inside the Problem structure.

[L Smith]

Thanks Walter.

I have coded the problem by defining a separate function to calculate the internal rate of return and defined it as a objective function in the script (see code below). I am unable to check the code because it makes a call to the function fcn2optimexpr which I believe exists only in the latest version of MATLAB. I will appreciate if some one can help with following questions:

1. Is the code below correct?
2. Is there any alternative formulation that can be used to avoid use of fcn2optimexpr?
3. Given that the objective is maximizing internal rate of return on the portfolio, does this mean the optimization problem is now non-linear? If so, should I be using fsolve instead of solve? Is there a way of making it linear to reduce run time?

clear

% Define data

prices = [99.74 91.22 98.71 103.75 97.15];

cashFlows = [4 5 2.5 5 4; 4 5 2.5 5 4; 4 5 2.5 5 4; 4 5 2.5 5 4; 4 5 102.5 5 4;4 5 0 105 104;

4 105 0 0 0; 104 0 0 0 0];

obligations = [5 7 7 6 8 7 20 0]';

%Create a cash flow matrix for internal rate of return calculation

cashFlows2 = [-prices; cashFlows]

nt=size(cashFlows,1)

nb=size(cashFlows,2)

bonds = optimvar('bonds',nb,'Type','integer','LowerBound',0,'UpperBound',1);

obj = fcn2optimexpr(@LiabDiscountR,bonds);

matchprob = optimproblem('Objective',obj);

matchprob.ObjectiveSense = 'maximize';

matchprob.Constraints.cover = cashFlows*bonds >= obligations;

showproblem(matchprob)

milpSol = solve(matchprob);

milpSol.bonds

function f = LiabDiscountR(bonds)

A = [cashFlows2*bonds]';

f = irr(A);

end

[Mary Fenelon]

Yes, it's nonlinear and a nonlinear problem with integer variables is much more difficult than a linear one.

First to get your code working. It's good practice to test your function before using it in fcn2optimexpr. If you did, you would get error messages that indicate you should

• Add cashFlows2 as an argument to LiabDiscountR
• Remove the transpose on A

Making those corrections and adding an initial point, you'll then get an error saying that nonlinear problems with integer variables aren't supported with OptimizationProblem.

One option is to solve with continuous bond variables, fix the bond variables that have integral or near-integral values by setting their bounds to that value, and re-solve with the remaining variables until all variable have integral values.

Another option is to solve it with the genetic algorithm, the ga function in Global Optimization Toolbox. It will be more successful if you provide it with an initial population of integral solutions. You could give it ones you obtain by variations on a fix-and-solve heuristic as above. You need to convert your problem to solver-based form to do this with the prob2struct and then set up the problem structure for what ga needs. I did that, just giving it one initial point:

%% 
% CashFlowMatchingIRR

clear
% Define data 
prices = [99.74 91.22 98.71 103.75 97.15];
cashFlows = [4 5 2.5 5 4; 4 5 2.5 5 4; 4 5 2.5 5 4; 4 5 2.5 5 4; 4 5 102.5 5 4;4 5 0 105 104;
4 105 0 0 0; 104 0 0 0 0];
obligations = [5 7 7 6 8 7 20 0]';
%Create a cash flow matrix for internal rate of return calculation 
cashFlows2 = [-prices; cashFlows]
nt=size(cashFlows,1)
nb=size(cashFlows,2)
bonds = optimvar('bonds',nb,'Type','integer','LowerBound',0,'UpperBound',1);

initialPt.bonds = ones(nb,1);
LiabDiscountR(initialPt.bonds,cashFlows2)
obj = fcn2optimexpr(@LiabDiscountR,bonds,cashFlows2);

matchprob = optimproblem('Objective',obj);
matchprob.ObjectiveSense = 'maximize';
matchprob.Constraints.cover = cashFlows*bonds >= obligations;
showproblem(matchprob)

gp = prob2struct(matchprob,initialPt);
gp.fitnessfcn = @(x)gp.objective(x(:));
gp.solver = 'ga';
gp.nvars = nb;
gp.options = optimoptions('ga');
gp.options.InitialPopulationMatrix = gp.x0';
[x,fval,exitflag,output] = ga(gp);

%%
function f = LiabDiscountR(bonds,cashFlows2)
A = cashFlows2*bonds;
f = irr(A);
end

There are some nonlinear integer solvers based on fmincon on File Exchange. There are also commercial and open source solvers for these problems and many have MATLAB interfaces.

[L Smith]

Thank you Mary. Presumably, this problem cannot be solved without fcn2optimexpr function?

[Mary Fenelon]

IRR is a nonlinear function--the only way to solve as a MILP is to approximate it with a linear function. One possibility is to take the IRR for each bond and use that as a linear objective coefficient. Is that a reasonable approximation or not? You'll need to decide that.

[L Smith]

Taking individual bonds irr to approximate portfolio irr could be one way of making the problem linear. In the code above that you kindly provided, it starts with the problem based approach and then converts to a Solver-based approach. As the objective function is not linear it uses fcn2optimexpr function. As far I can see this function is required if a problem based approach is used (as we start of in this case) and not for Solver-based. Wouldn’t it be better to only use the Solver based approach in this case from the start so that use of fcn2optimexpr will not be necessary?

[Mary Fenelon]

Sure, for this problem, going with Solver-based makes sense. There's only one group of variables, bonds, and it is one dimensional. There is only one group of constraints. You can use MATLAB matrix functions to write the constraints.

Problem-based is a win when there are several groups of variables, when the variables are multi-dimensional., or when there are several groups of constraints. For those sorts of problems, Problem-based is easier is easier for me to write the problem and easier for me to remember what I had in mind when I return to it after weeks or months.

[Walter Roberson]

Sometimes it is easier to write the expressions using the problem based approach, so that you can write in terms of variable names that are meaningful for you. However the problem based approach is not more powerful than the solver based approach: in most cases they use the same internal logic, sometimes simply converting to solver based internally and using that.

As you have found, the problem based approach is not yet as flexible as the solver based approach.

In setting up systems such as these, often the "better" approach is the one that makes it most likely that what is optimized is what you need to be optimized. The solver based approach is more flexible but it does require a lot of packing and unpacking of variables into a vector of values to be optimized over, and it can be a lot easier in such a way of writing that code to accidentally code values into the wrong place in arrays.

[L Smith]

Thank you Mary and Walter much appreciated. I agree with you that it is much easier to apply the problem based approach.