This example shows how to use the problem-based approach 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. For the solver-based approach, see Maximize Long-Term Investments Using Linear Programming: Solver-Based.

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 a fixed interest rate, compounds the investment 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 at the time of the investment.

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 using the optimization problem approach.

Start with a small example:

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 beginning of 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`MaturityYears`

represents the end of the year that bond`k`

is available for sale.Entry

`k`

of vector`InterestRates`

represents the percentage 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]; % Bond sale times MaturityYears = PurchaseYears + Maturity - 1; % Interest rates in percent InterestRates = [0;0;0;0;0;2;4;6;6]; % Return after one year of interest rt = 1 + InterestRates/100; plotInvestments(N,PurchaseYears,Maturity,InterestRates)

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(k)(1+{\rho}_{k}/100{)}^{{m}_{k}}.$$

Define $${\beta}_{k}=1+{\rho}_{k}/100$$ and define $${r}_{k}$$ as the total return of bond `k`

:

$${r}_{k}=(1+{\rho}_{k}/100{)}^{{m}_{k}}={\beta}_{k}^{{m}_{k}}.$$

x = optimvar('x',nPtotal,'LowerBound',0); % Total returns r = rt.^Maturity;

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}$$

Create an optimization problem for maximization, and include the objective function.

interestprob = optimproblem('ObjectiveSense','maximize'); interestprob.Objective = x(5)*r(5) + x(7)*r(7) + x(8)*r(8);

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}$$

investconstr = optimconstr(T,1); investconstr(1) = x(1) + x(6) == Capital_0; investconstr(2) = x(2) + x(8) + x(9) == r(1)*x(1); investconstr(3) = x(3) == r(2)*x(2); investconstr(4) = x(4) == r(3)*x(3); investconstr(5) = x(5) + x(7) == r(4)*x(4) + r(6)*x(6) + r(9)*x(9); interestprob.Constraints.investconstr = investconstr;

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 on the solution vector $$x$$. There is no explicit upper bound on the solution vector.

x.LowerBound = 0;

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'); [sol,fval,exitflag] = solve(interestprob,'options',options)

Solution found during presolve.

`sol = `*struct with fields:*
x: [9x1 double]

fval = 1.2625e+03

exitflag = OptimalSolution

The exit flag indicates that the solver found an optimal solution. The value `fval`

, returned as the second output argument, corresponds to the final wealth. Look at the final sum of investments, and the investment allocation 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,sol.x)

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\times {Capital}_{0}\\ {x}_{2}& \le Pmax\times ({\beta}_{1}{x}_{1}+{\beta}_{6}{x}_{6})\\ {x}_{3}& \le Pmax\times ({\beta}_{2}{x}_{2}+{\beta}_{6}^{2}{x}_{6}+{\beta}_{8}{x}_{8}+{\beta}_{9}{x}_{9})\\ {x}_{4}& \le Pmax\times ({\beta}_{3}{x}_{3}+{\beta}_{6}^{3}{x}_{6}+{\beta}_{8}^{2}{x}_{8}+{\beta}_{9}^{2}{x}_{9})\\ {x}_{5}& \le Pmax\times ({\beta}_{4}{x}_{4}+{\beta}_{6}^{4}{x}_{4}+{\beta}_{8}^{3}{x}_{8}+{\beta}_{9}^{3}{x}_{9})\\ {x}_{6}& \le Pmax\times {Capital}_{0}\\ {x}_{7}& \le Pmax\times ({\beta}_{4}{x}_{4}+{\beta}_{6}^{4}{x}_{4}+{\beta}_{8}^{3}{x}_{8}+{\beta}_{9}^{3}{x}_{9})\\ {x}_{8}& \le Pmax\times ({\beta}_{1}{x}_{1}+{\beta}_{6}{x}_{6})\\ {x}_{9}& \le Pmax\times ({\beta}_{1}{x}_{1}+{\beta}_{6}{x}_{6})\end{array}$$

```
% Maximum percentage to invest in any bond
Pmax = 0.6;
constrlimit = optimconstr(nPtotal,1);
constrlimit(1) = x(1) <= Pmax*Capital_0;
constrlimit(2) = x(2) <= Pmax*(rt(1)*x(1) + rt(6)*x(6));
constrlimit(3) = x(3) <= Pmax*(rt(2)*x(2) + rt(6)^2*x(6) + rt(8)*x(8) + rt(9)*x(9));
constrlimit(4) = x(4) <= Pmax*(rt(3)*x(3) + rt(6)^3*x(6) + rt(8)^2*x(8) + rt(9)^2*x(9));
constrlimit(5) = x(5) <= Pmax*(rt(4)*x(4) + rt(6)^4*x(6) + rt(8)^3*x(8) + rt(9)^3*x(9));
constrlimit(6) = x(6) <= Pmax*Capital_0;
constrlimit(7) = x(7) <= Pmax*(rt(4)*x(4) + rt(6)^4*x(6) + rt(8)^3*x(8) + rt(9)^3*x(9));
constrlimit(8) = x(8) <= Pmax*(rt(1)*x(1) + rt(6)*x(6));
constrlimit(9) = x(9) <= Pmax*(rt(1)*x(1) + rt(6)*x(6));
interestprob.Constraints.constrlimit = constrlimit;
```

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.

`[sol,fval] = solve(interestprob,'options',options);`

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 selected value of 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,sol.x)

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.

% 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; % Return after one year of interest rt = 1 + InterestRates/100; % Compute the return at the end of the maturity period for each bond: r = rt.^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 at the end of the year MaturityYears = PurchaseYears + Maturity - 1;

Compute the times when bonds can be bought or sold. The `buyindex`

matrix holds the potential purchase times, and the `sellindex`

matrix holds the potential sales times for each bond.

buyindex = false(nPtotal,T); % allocate nPtotal-by-T matrix for ii = 1:T buyindex(:,ii) = PurchaseYears == ii; end sellindex = false(nPtotal,T); for ii = 1:T sellindex(:,ii) = MaturityYears == ii; end

Set up the optimization variables corresponding to the bonds.

x = optimvar('x',nPtotal,1,'LowerBound',0);

Create the optimization problem and objective function.

interestprob = optimproblem('ObjectiveSense','maximize'); interestprob.Objective = sum(x(sellindex(:,T)).*r(sellindex(:,T)));

For convenience, create a temporary array xBuy, whose columns represent the bonds we can buy at each time period.

xBuy = repmat(x,1,T).*double(buyindex);

Similarly, create a temporary array xSell, whose columns represent the bonds we can sell at each time period.

xSell = repmat(x,1,T).*double(sellindex);

The return generated for selling these bounds is

xReturnFromSell = xSell.*repmat(r,1,T);

Create the constraint that the amount you invest in each time period is the amount that you sold in the previous time period.

interestprob.Constraints.InitialInvest = sum(xBuy(:,1)) == Capital_0; interestprob.Constraints.InvestConstraint = sum(xBuy(:,2:T),1) == sum(xReturnFromSell(:,1:T-1),1);

Solve the problem.

```
tic
[sol,fval,exitflag] = solve(interestprob,'options',options);
```

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 selected value of the constraint tolerance.

toc

Elapsed time is 0.196197 seconds.

How well did the investments do?

fprintf('After %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

To create constraints that limit the fraction of investments in each asset, set up a matrix that keeps track of the active bonds at each time. To express the constraint that each investment must be less than `Pmax`

times the total value, set up a matrix that keeps track of the value of each investment at each time. For this larger problem, set the maximum fraction that can be held to 0.4.

Pmax = 0.4;

Create an `active`

matrix corresponding to times when a bond can be held, and a `cactive`

matrix that holds the cumulative duration of each active bond. So the value of bond `j`

at time `t`

is `x(j)*(rt^cactive)`

.

active = double(buyindex | sellindex); for ii = 1:T active(:,ii) = double((ii >= PurchaseYears) & (ii <= MaturityYears)); end cactive = cumsum(active,2); cactive = cactive.*active;

Create the matrix whose entry (j,p) represents the value of bond j at time period p:

bondValue = repmat(x, 1, T).*active.*(rt.^(cactive));

Determine the total value of the investments at each time interval so you can impose the constraint on limited holdings. `mvalue`

is the money invested in all the bonds at the end of each time period, an `nPtotal`

-by-`T`

matrix.moneyavailable is the sum over the bonds of the money invested at the beginning of the time period, meaning the value of the portfolio at each time.

constrlimit = optimconstr(nPtotal,T); constrlimit(:,1) = xBuy(:,1) <= Pmax*Capital_0; constrlimit(:,2:T) = xBuy(:,2:T) <= repmat(Pmax*sum(bondValue(:,1:T-1),1), nPtotal, 1).*double(buyindex(:,2:T)); interestprob.Constraints.constrlimit = constrlimit;

Solve the problem with limited holdings.

```
tic
[sol,fval,exitflag] = solve(interestprob,'options',options);
```

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 selected value of the constraint tolerance.

toc

Elapsed time is 2.088101 seconds.

fprintf('After %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

To speed up the solver, try the dual-simplex algorithm.

options = optimoptions('linprog','Algorithm','dual-simplex'); tic [sol,fval,exitflag] = solve(interestprob,'options',options);

Optimal solution found.

toc

Elapsed time is 0.639190 seconds.

fprintf('After %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 less time to obtain the same solution.

To get a feel for the solution, 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,sol.x,false)