The Optimal Retirement Planner (ORP) contains a Monte Carlo linear programming (LP) implementation of a retirement planner

Abstract

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.

Introduction

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’ [3] last optimizer. Bosc was originally implemented to solve feed blending models for Precision Soya Inc in New Castle, Indiana.

Architecture

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 this

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

EXIT

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 the matrix.
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 follows:

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

EXIT

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:

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 matrix.

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 element:

Matrix(1,k) is the row index for element k.
Matrix(2,k) is the index into the coefficient pool for element k.

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 [2] 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:

-v 1

-v …

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 row name

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

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 ) = InvRet(x)*(AnnualRet(z)/0.08)

The justification for this is to be found in Appendix D.

MCDecide

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 Report.

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.

Conclusion

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 Investment Policy

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 high volatility.

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).

References:

Monte Carlo working paper; i-orp.com/ModelDescription/MonteCarlo.html.)
Denis Rarick, WHIZARD User Manual, Management Science Systems, Inc. 1977

###

A Monte Carlo Linear Programming Implementation

Abstract

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.

Introduction

Appendix C: Computational Algorithm

Appendix D: Incorporating the User's Investment Policy