Most Linear Programming (LP)
models include exogenous parameters whose values cannot be influenced by the models’
users. Commodity prices, interest rates and return on investments are examples
of exogenous parameters. Monte Carlo Linear Programming (MCLP) is a set of
instances for a particular model in which the values for the exogenous
parameters are generated randomly. The results of an MCLP is a statistical
report summarizing the behavior of the runs.
This paper presents the
architecture and implementation for extending a linear programming optimizer to
implement a Monte Carlo linear programming optimizer. An example of an
operational MCLP is given.
A conventional Linear Programming
(LP) model provides a single answer for a given set of parameters. Varying
those parameters for a set of solutions is called sensitivity analysis and
gives a picture of how the model behaves as the parameters change. This is a
very useful tool for the real world practitioner.
Some of the model’s
parameters are exogenous variables in that they are determined by environmental
factors outside of the practitioner’s control. Investment returns, interest
rates and crude oil prices are examples of exogenous variables. MCLP solves many
instances of the model, each with its exogenous variable randomly generated.
Whereas a conventional LP reports a solution view of the model’s variables, MCLP’s
output presents a statistical view of the behavior of the model over the range
of these randomly generated variables. Generating these random values is an
art form unto itself, peculiar to each application.
This paper presents a
generalized architecture for embedding the code for a MCLP into a conventional
LP system. Included is an example of a successfully application of this
technique in the form of an Internet based retirement decision support system,
the Optimal Retirement Planner (ORP, http://i-orp.com).
ORP is a linear
programming (LP) implementation of a retirement planner. It has been available,
at no charge,
on the Internet free since 1996. The purpose of ORP is demonstrate how a
computationally intensive process can be made accessible to the public.
The Monte Carlo option is a recent addition to ORP.
The MCLP is an extension to the
Bosc optimizer used to solve the ORP LP. Bosc is William Orchard-Hays’  last
optimizer. Bosc was originally implemented to solve feed blending models for Precision
Soya Inc in New Castle, Indiana.
Bosc is comprised of a main program and a library of
functions. The main program is available in source form so that it can be
modified. Several commercial LP systems are available in a form similar to
Abstractly the main program looks something like this:
CALL INPUT( … ) !Input
and set up the problem, from a standard MPS input file.
CALL SOLVE(…) !Get an
optimal or failed solution
CALL OUTPUT(…) !Generate
the standard MPS solution file
Embedding MCLP into Bosc requires adding 3 new functions to
the Bosc library and adjusting the main program. The three new functions are:
- MCBuild: analyze the LP matrix and set up the
infrastructure for inserting the random values into the matrix.
- MCFill: generate the random values and insert them into
- MCDecide: evaluate the most recent solution to decide
whether to go on for another Monte Carlo iteration or stop. If the
decision is to stop, write a file containing the summary of all runs.
These three functions are called from the main program as
shown below. This assumes that a control value came into the main program in
its command line. The settings for the control value are:
- Normal LP, no Monte Carlo processing.
- Debug mode, make one Monte Carlo iteration and stop.
- Monte Carlo mode, make multiple LP runs each with
different randomly generated values for certain exogenesis values.
The pseudo code logic for the revised main program is as
CALL INPUT(…) !Input and
set up the problem, from a standard MPS input file.
IF( control != 1 ) CALL MCBuild(…)
!Create Monte Carlo structures
100 CONTINUE !Loop
back here on each Monte Carlo iteration
IF( control != 1 ) CALL MCFill(…) !Generate
the random values
CALL SOLVE(…) !Get
an optimal or failed solution
IF( problem not optimal ) GOTO 100
!Generate new values and try again
IF( control = 3 and MCDecide( ) =0
) GOTO 100 !Not done yet
CALL OUTPUT(…) !Generate
the standard MPS solution file
When SOLVE is called the basis from the previous run is the
starting basis for the next solution. Only a few LP iterations are required to
reach the next optimal solution.
The three new functions assume that the in-memory LP matrix,
its descriptive structures and the optimizer’s solution vector are accessible.
All of the processing done by the functions is done from the optimizer’s
structures or from supporting data loaded by MCBuild. Beyond these changes to
the main program there are no other modifications to the Bosc LP system.
Monte Carlo Function Descriptions
This section describes the three functions that implement Monte
Carlo LP in Bosc. Each is described generally as a generic MCLP application
and then specifically as the ORP implementation.
MCBuild scans through the LP matrix, finds the coefficients
to be modified and sets up a structure for MCFill to use during Monte Carlo
iterations. MCBuild is called only once.
The LP matrix is described in accessible values and arrays:
m is the number of rows in the model.
n is the number of columns in the model
mn is the number of non null coefficients in the
RowNames (m) contains an 8 character name for each
row in the matrix.
ColNames(n) contains an 8 character name for each
column in the matrix.
ORP is a MathPro generated matrix. MathPro views an LP
matrix as row strips, that are collections of equations of identical form, and
column strips, that are collections of columns of identical structure. The
first character of each row name identifies its row strip with an alphabetic
character, A, B, C, …. Similarly for column strips.
ColDir(n+1) contains the integer index of the column’s
elements in the matrix array. ColDir(n+1) points the matrix cell after the last
matrix element. Thus the matrix entries for column j can be found in the elements
ColDir(j) through ColDir(j+1)-1. The number of elements in the column is
computed as ColDir(j+1) – ColDir(j).
Matrix(2,mn) contains two values for each matrix
- Matrix(1,k) is the row index for element k.
- Matrix(2,k) is the index into the coefficient pool for
Pool(nc) is a single precision, floating point array
of unique model coefficients. Most of the coefficients in ORP are 1, -1 or two
digit fractions. Bosc uses WHIZARD’s  coefficient pool implementation.
While there are nc coefficients in Pool the actual array is much bigger and it
had better be because MCBuild is going to need it.
MCLP has a few variables and arrays that are private and not
of interest to the LP side of things.
y is the number of years that are being modeled by
ORP. For a 55 year old retiree expecting to live until 90 y = 45, 10 years of
work left and 35 years in retirement.
InvRet(y) contains the investment returns specified
by the user in the ORP input form and put into the matrix by MathPro. These
values are used by MCFill while calculating the randomized return values
MCBuild is interested in columns in the column strips B, L
and Q as well as the rows in row strips C, F and M. The coefficients in blocks
(B, C), (L, F) and (Q, M) are to be randomized. ORP is a time dynamic model
and the structure of these blocks is:
The v’s are the coefficients that are the average return on
investment that are to be replaced by randomized values.
MCBuild iterates through ColName looking for names from the
target column strips. Columns for a strip are all grouped together which makes
things easier. When MCBuild finds the column strip of interest it scans through
the matrix elements, pointed to by ColDir, until it finds the row strip of
interest. The row index of the matrix element points to the row name in
RowNames. The rows of interest are identified by the first character of the
MCBuild moves the pool coefficient to InvRet(q). Then the
matrix pool index (Matrix(2,h) is changed to point to the next available pool
position Pool(nc+q). Then increment q, the count of new coefficients added to
the pool. Because of the structure of the block of coefficients, i.e. there is
only one coefficient of interest for each column; once a coefficient has been
processed we can move to the next matrix column.
y is the number of years being modeled. During the
processing of the last block only – (Q,M)- MCBuild counts the number of
intersections processed, which is also the number of years being modeled. This
count is of great interest to MCFill. ORP models the IRA and Roth IRA
separately for the retiree and spouse. The block (Q,M) represents the couples
after-tax account which is combined account. The column strips B and L contains
two matrix columns (retiree and spouse) for each year which is not too helpful.
InvRet, the last value of q and y are all kept in COMMON for
MCFill’s use. q = 5*y for a married couple and q = 3*y for a single retiree.
MCBuild finishes by initializing a hand full of variables
used by MCDecide. For example the count of Monte Carlo iterations is set to 0.
MCFill computes randomly generated coefficients and inserts
them into the matrix. MCFill has two steps. The first generates the random
returns for each year being modeled. The second scales that return by the
user’s volatility desires and stores it in the matrix pool.
AnnualRet(y) is local to MCFill. A random asset
return is generated for each year and stored in AnnualRet. The details of the
generation process can be found in Appendix B and Appendix C.
Next the InvRet array is scanned using index x running from
1 to q. The year being processed is:
z = mod(x-1, y)+1.
The new coefficient stored is
Pool(nc+x ) =
The justification for this is to be found in Appendix D.
MCDecide keeps statistics on the Monte Carlo process is
going and makes the decision as to when to stop.
S is amount of money, after taxes, that is available
for spending in the first year of retirement. (It is inflated annually
thereafter.) It is found in the last position of Bosc’s solution vector. (A back
up position is to scan the row names looking for its row index but since the
row is basic from the get go searching for it is usually not needed.) S is
ORP’s most important computed value, followed by the schedule in the Withdrawal
MCDecide keeps a running sum of S and the sum of the squares
of S. (The initial values of zero are set by MCBuild).
The average of S for all prior iterations is computed.
i, the Monte Carlo iteration count is incremented.
The sum and sum of the squares of S are incremented by S and S**2.
The new average is computed.
If the new average and the old average disagree by more than
the tolerance 0.1 and fewer than 500 iterations have been run then MCDecide
signals the main program to carry on and exits. 0.1 is actually $1,000 so
that’s not a very tight tolerance
Otherwise the process terminates. MCDecide writes a file in
the form of an html document containing the Monte Carlo statistics and exits. The
summary shows the number of Monte Carlo iterations, the mean value for S and
the first and second standard deviations below the mean and the worst trial
encountered. The html file is picked up by client and displayed to the user.
The process described here is specific to ORP. The
architecture is general and can be applied to any Monte Carlo LP
implementation. All that has to be done is to define the communicating common
and program the three functions.
Appendix B: Annual Return Domain
To create ORP's asset return domain the list of annual returns of the
S&P 500 stocks, sorted by year, is directed into two other lists.
List Pos contains returns where the previous year's return was positive.
i.e. Pos list contains r(j) if r(j-1) > 0.
Likewise the Neg list contains returns where the previous years return was
negative: i.e. Net list contains r(j) if r(j-1) <= 0.
These two lists are used by MCFill to create the randomized annual returns
as described in Appendix c.
Appendix C: Computational Algorithm
The idea is to compute a random annual investment return for each year j. ORP
does random sampling of the historical data stored in Pos and Neg. The
selected value is random, historically accurate and follows the trend set by
the previous coefficient.
Compute model's return m for year j as:
For j=1, randomly select m(1) from the combined list.
Set m(j) = Pos(k) if m(j-1) > 0 where k is a random integer 1<= k <=
501, j > 2.
Set m(j) = Neg(k) if m(j-1) < 0 where k is a random integer 1 <= k <.=
195, j > 2.
Appendix D: Incorporating the User's
From 1955 to 2009 the average annual return for the S&P 500 is 8%, or
.08. The user has specified her expected annual return on the ORP parameter
form is available in the LP matrix..
A specified return R < .08 is a conservative investment policy with a low
volatility. Expected returns R >= .08 are for an aggressive portfolio with a
Each m(j) is adjusted to reflect the users expectations (hopes). m(j) is
scaled to reflect this using the scale R/.08, where .08 is the average return
of the combine list of returns and R is the user specified expected asset returns.
Note that if R < .08 then R/.08 is less than 1.0 and return volatility is
reduced. Similarly for R >.08 the return volatility is increased. Therefore
the generated value is m`(j) = (R/.08) * m(j).
- Monte Carlo working paper; i-orp.com/ModelDescription/MonteCarlo.html.)
- Denis Rarick, WHIZARD User Manual, Management Science
Systems, Inc. 1977