# Maximize Long-Term Investments Using Linear Programming: Solver-Based

This example shows how to use the linprog solver in Optimization Toolbox® to solve an investment problem with deterministic returns over a fixed number of years T. The problem is to allocate your money over available investments to maximize your final wealth. This example uses the solver-based approach.

### Problem Formulation

Suppose that you have an initial amount of money Capital_0 to invest over a time period of T years in N zero-coupon bonds. Each bond pays an interest rate that compounds each year, and pays the principal plus compounded interest at the end of a maturity period. The objective is to maximize the total amount of money after T years.

You can include a constraint that no single investment is more than a certain fraction of your total capital.

This example shows the problem setup on a small case first, and then formulates the general case.

You can model this as a linear programming problem. Therefore, to optimize your wealth, formulate the problem for solution by the linprog solver.

• The starting amount to invest Capital_0 is $1000. • The time period T is 5 years. • The number of bonds N is 4. • To model uninvested money, have one option B0 available every year that has a maturity period of 1 year and a interest rate of 0%. • Bond 1, denoted by B1, can be purchased in year 1, has a maturity period of 4 years, and interest rate of 2%. • Bond 2, denoted by B2, can be purchased in year 5, has a maturity period of 1 year, and interest rate of 4%. • Bond 3, denoted by B3, can be purchased in year 2, has a maturity period of 4 years, and interest rate of 6%. • Bond 4, denoted by B4, can be purchased in year 2, has a maturity period of 3 years, and interest rate of 6%. By splitting up the first option B0 into 5 bonds with maturity period of 1 year and interest rate of 0%, this problem can be equivalently modeled as having a total of 9 available bonds, such that for k=1..9 • Entry k of vector PurchaseYears represents the year that bond k is available for purchase. • Entry k of vector Maturity represents the maturity period ${m}_{k}$ of bond k. • Entry k of vector InterestRates represents the interest rate ${\rho }_{k}$ of bond k. Visualize this problem by horizontal bars that represent the available purchase times and durations for each bond. % Time period in years T = 5; % Number of bonds N = 4; % Initial amount of money Capital_0 = 1000; % Total number of buying oportunities nPtotal = N+T; % Purchase times PurchaseYears = [1;2;3;4;5;1;5;2;2]; % Bond durations Maturity = [1;1;1;1;1;4;1;4;3]; % Interest rates InterestRates = [0;0;0;0;0;2;4;6;6]; plotInvestments(N,PurchaseYears,Maturity,InterestRates) ### Decision Variables Represent your decision variables by a vector x, where x(k) is the dollar amount of investment in bond k, for k=1..9. Upon maturity, the payout for investment x(k) is $x\left(k\right)\left(1+{\rho }_{k}/100{\right)}^{{m}_{k}}.$ Define ${\beta }_{k}=1+{\rho }_{k}/100$ and define ${r}_{k}$ as the total return of bond k: ${r}_{k}=\left(1+{\rho }_{k}/100{\right)}^{{m}_{k}}={\beta }_{k}^{{m}_{k}}.$ % Total returns finalReturns = (1+InterestRates/100).^Maturity; ### Objective Function The goal is to choose investments to maximize the amount of money collected at the end of year T. From the plot, you see that investments are collected at various intermediate years and reinvested. At the end of year T, the money returned from investments 5, 7, and 8 can be collected and represents your final wealth: $\underset{x}{\mathrm{max}}{x}_{5}{r}_{5}+{x}_{7}{r}_{7}+{x}_{8}{r}_{8}$ To place this problem into the form linprog solves, turn this maximization problem into a minimization problem using the negative of the coefficients of x(j): $\underset{x}{\mathrm{min}}{f}^{T}x$ with $f=\left[0;0;0;0;-{r}_{5};0;-{r}_{7};-{r}_{8};0\right]$ f = zeros(nPtotal,1); f([5,7,8]) = [-finalReturns(5),-finalReturns(7),-finalReturns(8)]; ### Linear Constraints: Invest No More Than You Have Every year, you have a certain amount of money available to purchase bonds. Starting with year 1, you can invest the initial capital in the purchase options ${x}_{1}$ and ${x}_{6}$, so: ${x}_{1}+{x}_{6}={Capital}_{0}$ Then for the following years, you collect the returns from maturing bonds, and reinvest them in new available bonds to obtain the system of equations: $\begin{array}{ll}{x}_{2}+{x}_{8}+{x}_{9}& ={r}_{1}{x}_{1}\\ {x}_{3}& ={r}_{2}{x}_{2}\\ {x}_{4}& ={r}_{3}{x}_{3}\\ {x}_{5}+{x}_{7}& ={r}_{4}{x}_{4}+{r}_{6}{x}_{6}+{r}_{9}{x}_{9}\end{array}$ Write these equations in the form $Aeqx=beq$, where each row of the $Aeq$ matrix corresponds to the equality that needs to be satisfied that year: $Aeq=\left[\begin{array}{ccccccccc}1& 0& 0& 0& 0& 1& 0& 0& 0\\ -{r}_{1}& 1& 0& 0& 0& 0& 0& 1& 1\\ 0& -{r}_{2}& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& -{r}_{3}& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& -{r}_{4}& 1& -{r}_{6}& 1& 0& -{r}_{9}\end{array}\right]$ $beq=\left[\begin{array}{c}{Capital}_{0}\\ 0\\ 0\\ 0\end{array}\right]$ Aeq = spalloc(N+1,nPtotal,15); Aeq(1,[1,6]) = 1; Aeq(2,[1,2,8,9]) = [-1,1,1,1]; Aeq(3,[2,3]) = [-1,1]; Aeq(4,[3,4]) = [-1,1]; Aeq(5,[4:7,9]) = [-finalReturns(4),1,-finalReturns(6),1,-finalReturns(9)]; beq = zeros(T,1); beq(1) = Capital_0; ### Bound Constraints: No Borrowing Because each amount invested must be positive, each entry in the solution vector $x$ must be positive. Include this constraint by setting a lower bound lb on the solution vector $x$. There is no explicit upper bound on the solution vector. Thus, set the upper bound ub to empty. lb = zeros(size(f)); ub = []; ### Solve the Problem Solve this problem with no constraints on the amount you can invest in a bond. The interior-point algorithm can be used to solve this type of linear programming problem. options = optimoptions('linprog','Algorithm','interior-point'); [xsol,fval,exitflag] = linprog(f,[],[],Aeq,beq,lb,ub,options); Solution found during presolve.  ### Visualize the Solution The exit flag is 1, indicating that the solver found a solution. The value -fval, returned as the second linprog output argument, corresponds to the final wealth. Plot your investments over time. fprintf('After %d years, the return for the initial$%g is $%g \n',... T,Capital_0,-fval); After 5 years, the return for the initial$1000 is $1262.48  plotInvestments(N,PurchaseYears,Maturity,InterestRates,xsol) ### Optimal Investment with Limited Holdings To diversify your investments, you can choose to limit the amount invested in any one bond to a certain percentage Pmax of the total capital that year (including the returns for bonds that are currently in their maturity period). You obtain the following system of inequalities: $\begin{array}{ll}{x}_{1}& \le Pmax×{Capital}_{0}\\ {x}_{2}& \le Pmax×\left({\beta }_{1}{x}_{1}+{\beta }_{6}{x}_{6}\right)\\ {x}_{3}& \le Pmax×\left({\beta }_{2}{x}_{2}+{\beta }_{6}^{2}{x}_{6}+{\beta }_{8}{x}_{8}+{\beta }_{9}{x}_{9}\right)\\ {x}_{4}& \le Pmax×\left({\beta }_{3}{x}_{3}+{\beta }_{6}^{3}{x}_{6}+{\beta }_{8}^{2}{x}_{8}+{\beta }_{9}^{2}{x}_{9}\right)\\ {x}_{5}& \le Pmax×\left({\beta }_{4}{x}_{4}+{\beta }_{6}^{4}{x}_{4}+{\beta }_{8}^{3}{x}_{8}+{\beta }_{9}^{3}{x}_{9}\right)\\ {x}_{6}& \le Pmax×{Capital}_{0}\\ {x}_{7}& \le Pmax×\left({\beta }_{4}{x}_{4}+{\beta }_{6}^{4}{x}_{4}+{\beta }_{8}^{3}{x}_{8}+{\beta }_{9}^{3}{x}_{9}\right)\\ {x}_{8}& \le Pmax×\left({\beta }_{1}{x}_{1}+{\beta }_{6}{x}_{6}\right)\\ {x}_{9}& \le Pmax×\left({\beta }_{1}{x}_{1}+{\beta }_{6}{x}_{6}\right)\end{array}$ Place these inequalities in the matrix form Ax <= b. To set up the system of inequalities, first generate a matrix yearlyReturns that contains the return for the bond indexed by i at year j in row i and column j. Represent this system as a sparse matrix. % Maximum percentage to invest in any bond Pmax = 0.6; % Build the return for each bond over the maturity period as a sparse % matrix cumMaturity = [0;cumsum(Maturity)]; xr = zeros(cumMaturity(end-1),1); yr = zeros(cumMaturity(end-1),1); cr = zeros(cumMaturity(end-1),1); for i = 1:nPtotal mi = Maturity(i); % maturity of bond i pi = PurchaseYears(i); % purchase year of bond i idx = cumMaturity(i)+1:cumMaturity(i+1); % index into xr, yr and cr xr(idx) = i; % bond index yr(idx) = pi+1:pi+mi; % maturing years cr(idx) = (1+InterestRates(i)/100).^(1:mi); % returns over the maturity period end yearlyReturns = sparse(xr,yr,cr,nPtotal,T+1); % Build the system of inequality constraints A = -Pmax*yearlyReturns(:,PurchaseYears)'+ speye(nPtotal); % Left-hand side b = zeros(nPtotal,1); b(PurchaseYears == 1) = Pmax*Capital_0; Solve the problem by investing no more than 60% in any one asset. Plot the resulting purchases. Notice that your final wealth is less than the investment without this constraint. [xsol,fval,exitflag] = linprog(f,A,b,Aeq,beq,lb,ub,options); Minimum found that satisfies the constraints. Optimization completed because the objective function is non-decreasing in feasible directions, to within the function tolerance, and constraints are satisfied to within the constraint tolerance.  fprintf('After %d years, the return for the initial$%g is $%g \n',... T,Capital_0,-fval); After 5 years, the return for the initial$1000 is $1207.78  plotInvestments(N,PurchaseYears,Maturity,InterestRates,xsol) ### Model of Arbitrary Size Create a model for a general version of the problem. Illustrate it using T = 30 years and 400 randomly generated bonds with interest rates from 1 to 6%. This setup results in a linear programming problem with 430 decision variables. The system of equality constraints is represented by a sparse matrix Aeq of dimension 30-by-430 and the system of inequality constraints is represented by a sparse matrix A of dimension 430-by-430. % for reproducibility rng default % Initial amount of money Capital_0 = 1000; % Time period in years T = 30; % Number of bonds N = 400; % Total number of buying oportunities nPtotal = N+T; % Generate random maturity durations Maturity = randi([1 T-1],nPtotal,1); % Bond 1 has a maturity period of 1 year Maturity(1:T) = 1; % Generate random yearly interest rate for each bond InterestRates = randi(6,nPtotal,1); % Bond 1 has an interest rate of 0 (not invested) InterestRates(1:T) = 0; % Compute the return at the end of the maturity period for each bond: finalReturns = (1+InterestRates/100).^Maturity; % Generate random purchase years for each option PurchaseYears = zeros(nPtotal,1); % Bond 1 is available for purchase every year PurchaseYears(1:T)=1:T; for i=1:N % Generate a random year for the bond to mature before the end of % the T year period PurchaseYears(i+T) = randi([1 T-Maturity(i+T)+1]); end % Compute the years where each bond reaches maturity SaleYears = PurchaseYears + Maturity; % Initialize f to 0 f = zeros(nPtotal,1); % Indices of the sale oportunities at the end of year T SalesTidx = SaleYears==T+1; % Expected return for the sale oportunities at the end of year T ReturnsT = finalReturns(SalesTidx); % Objective function f(SalesTidx) = -ReturnsT; % Generate the system of equality constraints. % For each purchase option, put a coefficient of 1 in the row corresponding % to the year for the purchase option and the column corresponding to the % index of the purchase oportunity xeq1 = PurchaseYears; yeq1 = (1:nPtotal)'; ceq1 = ones(nPtotal,1); % For each sale option, put -\rho_k, where \rho_k is the interest rate for the % associated bond that is being sold, in the row corresponding to the % year for the sale option and the column corresponding to the purchase % oportunity xeq2 = SaleYears(~SalesTidx); yeq2 = find(~SalesTidx); ceq2 = -finalReturns(~SalesTidx); % Generate the sparse equality matrix Aeq = sparse([xeq1; xeq2], [yeq1; yeq2], [ceq1; ceq2], T, nPtotal); % Generate the right-hand side beq = zeros(T,1); beq(1) = Capital_0; % Build the system of inequality constraints % Maximum percentage to invest in any bond Pmax = 0.4; % Build the returns for each bond over the maturity period cumMaturity = [0;cumsum(Maturity)]; xr = zeros(cumMaturity(end-1),1); yr = zeros(cumMaturity(end-1),1); cr = zeros(cumMaturity(end-1),1); for i = 1:nPtotal mi = Maturity(i); % maturity of bond i pi = PurchaseYears(i); % purchase year of bond i idx = cumMaturity(i)+1:cumMaturity(i+1); % index into xr, yr and cr xr(idx) = i; % bond index yr(idx) = pi+1:pi+mi; % maturing years cr(idx) = (1+InterestRates(i)/100).^(1:mi); % returns over the maturity period end yearlyReturns = sparse(xr,yr,cr,nPtotal,T+1); % Build the system of inequality constraints A = -Pmax*yearlyReturns(:,PurchaseYears)'+ speye(nPtotal); % Left-hand side b = zeros(nPtotal,1); b(PurchaseYears==1) = Pmax*Capital_0; % Add the lower-bound constraints to the problem. lb = zeros(nPtotal,1); ### Solution with No Holding Limit First, solve the linear programming problem without inequality constraints using the interior-point algorithm. % Solve the problem without inequality constraints options = optimoptions('linprog','Algorithm','interior-point'); tic [xsol,fval,exitflag] = linprog(f,[],[],Aeq,beq,lb,[],options); Minimum found that satisfies the constraints. Optimization completed because the objective function is non-decreasing in feasible directions, to within the function tolerance, and constraints are satisfied to within the constraint tolerance.  toc Elapsed time is 0.043860 seconds.  fprintf('\nAfter %d years, the return for the initial$%g is $%g \n',... T,Capital_0,-fval); After 30 years, the return for the initial$1000 is $5167.58  ### Solution with Limited Holdings Now, solve the problem with the inequality constraints. % Solve the problem with inequality constraints options = optimoptions('linprog','Algorithm','interior-point'); tic [xsol,fval,exitflag] = linprog(f,A,b,Aeq,beq,lb,[],options); Minimum found that satisfies the constraints. Optimization completed because the objective function is non-decreasing in feasible directions, to within the function tolerance, and constraints are satisfied to within the constraint tolerance.  toc Elapsed time is 1.055830 seconds.  fprintf('\nAfter %d years, the return for the initial$%g is $%g \n',... T,Capital_0,-fval); After 30 years, the return for the initial$1000 is $5095.26  Even though the number of constraints increased by an order of 10, the time for the solver to find a solution increased by an order of 100. This performance discrepancy is partially caused by dense columns in the inequality system shown in matrix A. These columns correspond to bonds with a long maturity period, as shown in the following graph. % Number of nonzero elements per column nnzCol = sum(spones(A)); % Plot the maturity length vs. the number of nonzero elements for each bond figure; plot(Maturity,nnzCol,'o'); xlabel('Maturity period of bond k') ylabel('Number of nonzero in column k of A') Dense columns in the constraints lead to dense blocks in the solver's internal matrices, yielding a loss of efficiency of its sparse methods. To speed up the solver, try the dual-simplex algorithm, which is less sensitive to column density. % Solve the problem with inequality constraints using dual simplex options = optimoptions('linprog','Algorithm','dual-simplex'); tic [xsol,fval,exitflag] = linprog(f,A,b,Aeq,beq,lb,[],options); Optimal solution found.  toc Elapsed time is 0.242389 seconds.  fprintf('\nAfter %d years, the return for the initial$%g is $%g \n',... T,Capital_0,-fval); After 30 years, the return for the initial$1000 is $5095.26  In this case, the dual-simplex algorithm took much less time to obtain the same solution. ### Qualitative Result Analysis To get a feel for the solution found by linprog, compare it to the amount fmax that you would get if you could invest all of your starting money in one bond with a 6% interest rate (the maximum interest rate) over the full 30 year period. You can also compute the equivalent interest rate corresponding to your final wealth. % Maximum amount fmax = Capital_0*(1+6/100)^T; % Ratio (in percent) rat = -fval/fmax*100; % Equivalent interest rate (in percent) rsol = ((-fval/Capital_0)^(1/T)-1)*100; fprintf(['The amount collected is %g%% of the maximum amount$%g '... 'that you would obtain from investing in one bond.\n'... 'Your final wealth corresponds to a %g%% interest rate over the %d year '... 'period.\n'], rat, fmax, rsol, T)
The amount collected is 88.7137% of the maximum amount \$5743.49 that you would obtain from investing in one bond. Your final wealth corresponds to a 5.57771% interest rate over the 30 year period. 
plotInvestments(N,PurchaseYears,Maturity,InterestRates,xsol,false) 