October 17, 2014

©
Minimizing the Probability of Ruin in Retirement

CHRISTOPHER J. ROOK*

ABSTRACT

Retirees who exhaust their savings while still alive are said to experience financial ruin. The
savings are typically grown during life’s accumulation phase then spent during the retirement
decumulation phase. Extensive research into invest-and-harvest decumulation strategies has been
conducted, but recommendations differ markedly. This has likely been a source of concern and
confusion for the retiree. Our goal is to find what has heretofore been elusive, namely an optimal
decumulation strategy. Optimality implies that no alternate strategy exists or can be constructed
that delivers a lower probability of ruin, given a fixed inflation-adjusted withdrawal rate.

* The author is a consultant statistical programmer and studies in the Department of Systems Engineering
at Stevens Institute of Technology. All appendices referred to in the text can be found in the Internet
Appendix document that accompanies this research.

The retirement literature comprises a vast body of research that has evolved disparately with
financial economists optimizing portfolio-choice or life-cycle models and practitioners adhering
to rules-of-thumb (Milevsky and Huang (2010)). Economists attempt to smooth a person’s
utility of consumption over their earning and retirement years and have criticized practitioner
heuristics as being ad-hoc, simplistic, inefficient, misguided, far from optimal, and potentially
unsafe (see Kotlikoff (2008); Scott, Sharpe, and Watson (2009); Irlam and Tomlinson (2014);
Finke, Pfau, and Blanchett (2013)). Practitioners counter that they consider risk tolerance when
forming and assessing a client’s financial plan (Mitchell and Smetters (2013)), and assert that
most clients seek constant, not variable spending and consumption during retirement. Another
practitioner concern is that utility functions are difficult to estimate accurately and if
misspecified results in maximization of the wrong expression (Kitces (2012)).
Our focus here is on the most ubiquitous of practitioner heuristics, namely the safe
withdrawal rate (WR). This heuristic posits that the retiree may withdraw WR-% from their
portfolio during year one and adjust it for inflation each year thereafter. The popularity of this
heuristic among practitioners is irrefutable. Scott et al. note that “retail brokerage firms, mutual
fund companies, retirement groups, investor groups, financial websites, and the popular financial
press all recommend it.” Milevsky and Huang (2010) add that “it is hard to overstate the
influence” this heuristic has had on modern retirement planning, noting that its results have been
reported “thousands of times in the last two decades” making it “destined for the same
immortality enjoyed by other rules-of-thumb”.
The effectiveness of a safe WR strategy is measured by its success rate which is the
percentage of times the portfolio lasts a fixed number of years, or outlives the retiree. Letting
P(Ruin) denote the probability of ruin, the success rate is 100%*(1-P(Ruin)). Maximizing the
success rate and minimizing P(Ruin) are thus equivalent optimization problems. Retirees using a
safe WR strategy face the difficult decision of setting its initial value and determining an asset
allocation that is appropriate.

Withdraw too little and face a reduced standard of living;

withdraw too much and face financial ruin. Likewise, investing the wrong proportion in stocks
or bonds can lead to ruin from excessive or insufficient risk. Unfortunately there is no consensus

1

regarding a precise “safe” initial WR or the optimal asset allocation for retirees. Generally a WR
between 3% and 6% is suggested but success rates vary considerably. The term glide-path refers
to asset allocations that either shed or acquire equities systematically over time. Some advocate
for a declining, and others for a rising glide-path during retirement. With such a large cohort (>
76 million U.S. Baby Boomers) on the precipice of spending their savings, finding optimal
decumulation strategies is increasingly important. The first Baby Boomers to retire at age 65
began doing so in 2011, and these retirements will continue at a rate of 10,000/day until 2030.1
The general approach to retirement planning by practitioners has been to develop
intuition-based heuristics and test them using historical data, bootstrapping, or simulation. These
heuristics usually apply to withdrawals or to the glide-path. Our goal here is to introduce a
theoretical framework for the safe WR heuristic into the literature. In addition, we optimize the
glide-path to achieve a minimum P(Ruin) for retirees. We show that, when used optimally, the
safe WR is not merely a heuristic but an academically sound approach to retirement planning
based on theory spanning a wide array of disciplines. At each time point during decumulation
there exists a minimum P(Ruin) based on an optimal asset allocation and we seek to find these
values. We formulate a model that reveals precisely these measures as a function of both time
and a quantity we term the ruin factor. Our model is optimal in the sense that no alternate WR
strategy can deliver a lower P(Ruin), but it is also intractable under common market
assumptions. We use numerical techniques to approximate the optimal solution and the user can
achieve any desired degree of accuracy by adjusting the discretization precision.
This research is organized as follows. In Section I the literature of practitioner-based
decumulation strategies is reviewed and in Section II we formulate and optimize models for
retirement horizons that are fixed or random in length. In Section III we implement these models
using historical market returns to estimate randomness, and we extend the random horizon model
to a group of retirees. In Section IV we conclude with a non-technical summary and provide
source code from a full implementation in the Appendix. Throughout the discussion we offer
suggestions for extending these models in various dimensions, for example by adding asset
classes, tailoring the expense ratio, and using forecasts to vary the return distribution across time.

2

I. Literature Review
Bengen (1994) determined the safe initial WR for retirees by analyzing U.S. historical
return sequences on portfolios of stocks and bonds. He concluded that a 4.15% inflationadjusted WR would not have been exhausted over any prior 30-year period using 50% stocks, and
therefore considered it safe. Bengen (1994) recommended a fixed equity allocation between
50% and 75% which would only change if the retiree’s objectives changed. This research led to
what is commonly referred to as the “4% rule”. Bengen (2006) updated this research to include
small cap stocks in the portfolio and used the term SAFEMAX to refer to the highest WR that
yields a >$0 terminal balance after 30 years for all historical sequences analyzed. When at least
20% of the portfolio’s stock was invested in small caps, the SAFEMAX grew to 4.42%.
Cooley, Hubbard, and Walz (1998) also applied historical return sequences to retirement
portfolios of stocks and bonds. They concluded that a 4% WR on accounts with 50% and 75%
stocks had success rates (>$0 terminal balance) of 95% and 98%, respectively over 30-year
horizons. In 2011, Cooley et al. updated the analysis using 14 years of new data and found that
portfolios with 50% and 75% stocks had success rates of 96% and >99%, respectively.
Pye (2000) analyzed the sustainable WR by making distributional assumptions about
investment returns and tested various withdrawal strategies using simulation. Portfolios were
constructed using stocks and TIPS and withdrawals were reduced once thresholds became
violated. If the portfolio approaches ruin withdrawals are restricted, reducing the probability of a
shortfall. Instead of measuring the probability of ruin, Pye (2000) measured the probability of
hitting the limits and facing a reduced standard of living (success=no reduction). A 4% WR
using 100% stocks and a 4.5% WR using 40% stocks both had 35-year success rates of 81%.
Guyton (2004) developed decision rules for making withdrawals and tested them on the
30-year period between 1973 and 2003, considered challenging due to high inflation and two
severe market downturns. Fixed portfolios having either 65% or 80% equities were used and the
rules limited inflation adjustments on WR to years when the account balance did not decline.
Withdrawals were also taken from stocks and bonds in a systematic manner that drew first from
accounts with outsized prior year gains. Using these rules, Guyton (2004) found that portfolio

3

equity ratios of 65% and 80% could sustain an initial WR of 5.4% and 5.8%, respectively for the
30-year period studied plus 10 additional years. The foregoing of an inflation-adjustment at
predetermined times thus results in a higher WR when compared to most other research.
Milevsky and Robinson (2005) introduced the stochastic present value (SPV) to derive
the ruin PDF over a retiree’s life expectancy which was approximated by an exponential RV.
The SPV is the present value of the retiree’s future spending discounted using random market
returns. Ruin occurs if the SPV exceeds the retiree’s starting account balance. Milevsky and
Robinson (2005) found that under the assumption of lognormal portfolio returns the SPV RV
follows a reciprocal gamma distribution with parameters directly related to the parameters of the
lognormal returns (mean and variance) and the exponential life expectancy (rate). The solution
is closed form in nature, but serves as an approximation to the true probability of ruin when
exponential lifetimes are used. Milevsky and Robinson (2005) found that equity ratios of 50%
and 100% yield success rates of 91.0% and 87.7%, respectively using a 4% WR.
Blanchett (2007) analyzed the safe WR by bootstrap sampling historical returns on
portfolios of stocks and bonds. Fixing the success rate at 95%, Blanchett (2007) found that a
portfolio using 40% equities had the highest WR of 4.2% over a 30-year horizon. Blanchett
(2007) also analyzed 4 downward sloping (linear, stair, concave, convex) equity glide-paths
along with the constant allocation (no glide-path). Each was tested with different starting equity
percentages and resulted in 43 separate paths analyzed on horizons of 20 to 40 years. Blanchett
(2007) found that a constant allocation performed best. Specifically, over 30-years a 4% WR
with constant allocation of 30% stocks would generate the lowest failure rate of 2.06%.
Stout (2008) formulated retirement ruin as the objective of a stochastic optimization
problem and solved it using simulation. Fixed allocation percentages to stocks and bonds were
varied until the minimum probability of ruin was found. The objective function thus depends on
random quantities such as market returns and age of death using survival probabilities from the
U.S. Centers for Disease Control. Stout (2008) found that a 65-year old retiree using a 4.25%
WR would achieve a success rate of 96.75%, assuming an optimal allocation of 50.25% in stocks.
Further, as WR increases the equity ratio must also increase to maintain optimality.

4

Target-date (TD) funds are financial products that implement declining equity glide-paths
based on the retirement date. Some cease at retirement, while others continue reallocating
through retirement. Morningstar, Inc. (2009) surveyed financial firms and summarized the glidepaths across a myriad TD funds. They found that retirees from 1995, 2000, 2010, and 2015
currently have on average 20%, 34.5%, 50.2%, and 62.1% in equities, respectively. This survey
provides a great deal of insight into the glide-paths used by financial firms in their TD funds. On
average the retiree begins with about 62% equities which decreases to 20% after two decades.
There is not widespread agreement that declining equity glide-paths outperform fixed
allocations, however. Similar to Blanchett (2007), Cohen, Gardner, and Fan (2010) found that a
constant equity mix will outperform a declining equity glide-path almost universally, raising
doubts about TD fund strategies. Cohen et al. (2010) used a fixed retirement horizon and found
that retirees should have fewer equities than is often recommended. Namely, equity percentages
of 32%, 60%, and 80% yielded success rates of 94%, 88%, and 84%, respectively using a 4%
WR. Their rational for a lower equity ratio is to prevent against sequence of return (SOR) risk,
which is the risk of losses early in retirement.

When coupled with inflation-adjusted

withdrawals, early losses increase the probability of experiencing financial ruin.
Fan, Murray, and Pittman (2013) developed an adaptive allocation model that tracks the
retiree’s account balance over time and adjusts the asset allocation to maximize expected wealth
over 20 years, prior to buying an annuity. A decision tree is built using combinations of market
returns and their impact on different asset allocations between time points. By folding back the
tree an optimal policy reveals itself and by eliminating sub-optimal nodes it serves as a guide for
the retiree over time. The allocation is optimized at each node in the tree using an objective
function that punishes shortfalls quadratically. The adaptive model results in higher expected
surpluses and lower expected losses compared to fixed strategies. A feature of this model is that
it becomes most conservative when current assets match discounted future spending. When the
retiree experiences strong returns the shortfall penalty approaches zero revealing a feedback
mechanism that encourages the building of more wealth as wealth grows. Therefore the adaptive
model advocates for a higher equity ratio in both underperforming and outperforming markets.

5

Pfau and Kitces (2014) researched various equity glide-paths for periods of typical and
atypical (sub-standard) returns using simulation. For 2 of the 3 periods tested, a 4% WR yielded
the highest success rate using equity glide-paths that start between 20% and 40%, and end
between 60% and 80%. These glide-paths also led to the smallest shortfall magnitude for all
periods analyzed. Rising equity glide-paths thus outperformed both static and declining glidepaths. Pfau and Kitces (2014) offer the intuition behind these findings, namely that the portfolio
assumes greatest SOR risk at the beginning when it is largest, and the retiree should limit their
equity exposure. Over time this risk lessens and adding equities is warranted. Pfau and Kitces
(2014) note that a rising glide-path during a poor market allows the retiree to accumulate stocks
when prices are low, reaping the benefits of a market upturn. Accumulating equities during a
rising market is viewed as an opportunity to build a legacy for heirs.
II. Formulating and Solving the Model
We formulate an optimal retirement decumulation model in Section II and implement it
in Section III. The solution can serve as a roadmap for decumulating retirees. Most conclusions
are justified in appendices and the reader is encouraged to review those sections as well.
A. Definitions and Notation
Assume the retiree owns an investment account consisting of stocks and bonds and that
financial ruin occurs iff the retiree empties their account while still alive. Define:
· t = time point (t=0, 1, …, TD), with retirement starting at t=0
· TD = time of last withdrawal prior to death (can be either fixed or random)
· TD-k = k time points before TD after making the withdrawal (there are k
withdrawals remaining at time t=TD-k)
· Ruin(t) = the event of ruin is experienced at time t
· RuinC(t) = the event of ruin is not experienced at time t
· $A = value of retirement account at time t=0
· α = proportion of retirement account in stocks (0 ≤ α ≤ 1)
· (1-α) = proportion of retirement account in bonds

6

· WR = initial withdrawal rate adjusted for inflation at each time point
· It = inflation rate between times t-1 and t
· R(t,α) = total rate of return (RoR) for the retirement account with 100*α% stocks and
100*(1-α)% bonds between times t-1 and t
· r(t,α) | It, R(t,α) = inflation-adjusted RoR for the retirement account with 100*α%
stocks and 100*(1-α)% bonds between times t-1 and t, given It and R(t,α)
· ER = expense ratio charged by the financial institution per time t (fixed or α*ER(s) +
(1-α)*ER(b) where ER(s), ER(b) are stock and bond expense ratios)
· ȓ(t,α) = (1-ER)*(1+r(t,α) | It, R(t,α)) inflation/expense-adjusted compounded return on
the account with 100*α% stocks and 100*(1-α)% bonds between times t-1 and t
The quantities It, R(t,α), r(t,α) | It, R(t,α), and ȓ(t,α) are continuous RVs and TD can be fixed or
(discrete) random. WR is based on $A and is withdrawn at time t after adjusting for inflation.
The real RoR for an α-based portfolio at time t is r(t,α) | It, R(t,α), which is calculated after R(t,α) and
It are observed, therefore it is defined conditionally. The quantity ȓ(t,α) reflects the inflation/
expense-adjusted compounded return for an α-based portfolio at time t. The RV ȓ(t,α) thus
depends on 5 quantities: t, It, R(t,α), α, ER, where t, α, and ER are deterministic and It, R(t,α)
random. The real balance at time t-1 multiplied by ȓ(t,α) is the real balance at time t. Finally,
Ruin(t) and RuinC(t) are random events that occur with some probability.
The event Ruin(t) occurs iff the account balance cannot sustain the withdrawal at time t.
It can only occur once. When expressing P(Ruin(t)) we condition on no ruin occurring before
time t. For example, the probability of ruin at time t=3 is P(Ruin(3)| RuinC(1) ∩RuinC(2)). The
conditioning component can be expressed as P(RuinC(1) ∩ RuinC(2)) = P(RuinC(1))*P(RuinC(2) |
RuinC(1)). The probability of experiencing Ruin(3) is then expressed as:
∩

P Ruin 3 | Ruin 1 ∩ Ruin 2

∗

∩

(1)

The probability of ruin, P(Ruin), changes over time. If the retiree using a fixed WR owns
an α-based portfolio whose real account balance doubles between times t=0 and t=1 then a very
different P(Ruin) exists at time t=1, than did at time t=0. We assume the retiree’s objective is to

7

minimize P(Ruin), not build wealth, given WR. Ruin probabilities can thus perish quickly and
should be reevaluated at regular intervals by the retiree.
B. The Retirement Sequence
We assume the retiree starts time t=0 with an account balance of $A. The withdrawal
rate WR is based on $A and adjusted for inflation at each time point. At time t=1 the retiree
attempts their first withdrawal of ($A)*(WR)*(1+I1) from the account which has compounded
returns for one time point.

If the account cannot sustain the withdrawal Ruin(1) occurs,

C

otherwise Ruin (1) occurs. The pre-withdrawal account balance is ($A)*(1+R(1,α))*(1-ER), where
R(1,α) is the portfolio’s total RoR between times t=0 and t=1.

Under an optimal strategy,

assuming RuinC(1), the retiree should reassess P(Ruin) at any future time point (t=2, 3, …, TD)
and reallocate the portfolio’s assets optimally.
This scenario assumes no withdrawal is made at time t=0. We base the model on this
assumption and refer to it as standard form. In standard form, the first asset allocation decision
occurs at time t=0 and the first withdrawal occurs at time t=1. If a retiree presents differently
they will be converted to standard form before applying the model. For example, a retiree with
time t=0 balance of $B that requires funding between times t=0 and t=1 can be converted to
standard form by using a time t=0 withdrawal rate of W0=WR/(1+WR) applied to $B. The
withdrawal at time t=0 is (W0)*($B) and the remaining balance is $A=$B*(1-W0). The standard
form model is then applied using $A and WR. Note that this conversion preserves the retiree’s
purchasing power from time t=0 to time t=1 since (W0)*($B)*(1+I1) = (WR)*($A)*(1+I1).2
C. The Ruin Factor
We now define a quantity that encapsulates much information needed to assess P(Ruin).
It should be updated by the retiree at each time point and is actionable in the sense that it helps
reveal the precise asset allocation that minimizes P(Ruin) for the remainder of retirement. Let
r(t,α) denote r(t,α) | It, R(t,α). Define the ruin factor, RF(t) as WR for t = 0 and:
RF(t) = RF(t-1) / [(1+r(t,α))(1-ER) – RF(t-1)]
= RF(t-1) / [ȓ(t,α) – RF(t-1)], for t = 1, 2, …, TD.

8

(2a)
(2b)

Given no ruin prior to time t, we show that the event Ruin(t) occurs as:
Ruin(t) | RuinC(1) ∩ … ∩ RuinC(t-1) ↔ (1+ r(t,α))*(1-ER) ≤ RF(t-1)
↔ ȓ(t,α) ≤ RF(t-1)

(3a)
(3b)

(See Appendix B.) At time t=t-1, RF(t-1) is known and Ruin(t) occurs iff ȓ(t,α) ≤ RF(t-1).
Assuming the CDF of ȓ(t,α) is known or estimated, RF(t-1) can be calculated and P(Ruin(t))
evaluated across various α at time t=t-1. RF(t) can be defined using any unit of time assuming
ȓ(t,α) is the appropriate return. Note that RF(t) is positive until ruin occurs (if it occurs) when it
becomes negative (or ∞). To see this note that at any time t, RuinC(t) occurs ↔:
ȓ(t,α) > RF(t-1)

(4)

This condition states that the upcoming time point’s return ȓ(t,α), must exceed RF(t-1) to avoid
ruin. If this fails to hold, ȓ(t,α) ≤ RF(t-1) and RF(t) is either negative (when ȓ(t,α) < RF(t-1)) or
undefined (when ȓ(t,α) = RF(t-1)). If RF(t) declines between time points the retiree has improved
their position with respect to ruin, and vice versa. RF(t) declines when the denominator exceeds
1, namely ȓ(t,α) > 1 + RF(t-1). If the retiree generates returns (r(t,α) - ER - r(t,α)*ER) that exactly
equal WR=RF(0), then RF(t) is constant and equals WR over time.
If RuinC(t) occurs the real account balance at time t is ($A)*RF(0)/RF(t). (See Appendix
A.) This reveals that the reciprocal of RF(t) equals the number of real withdrawals remaining at
time t. At time t, RF(t) is computed and known. Future ruin factors are unknown since they are
functions of unobserved market returns, thus they are RVs. At the start of retirement RF(0)=WR
is the only known ruin factor. The retiree exercises some control over future RF(t) values via the
choice of α since RF(t)=RF(t-1)/[ȓ(t,α) – RF(t-1)]. If a low-α portfolio is used between times t-1
and t, the retiree is more confident of ȓ(t,α), and thus RF(t). Conversely, if a high-α portfolio is
used ȓ(t,α) has greater uncertainty and so too does RF(t).
D. Example of Two Retirees Facing Different SORs
Consider two retirees with rebalanced α=0.5 portfolios. Let the unit of time be years with
TD=30. Both begin time t=0 with account balances of $A use WR=4% with ER=0.5%. To avoid
ruin, the retirees must have >$0 account balances after making the time TD=30 withdrawal. At

9

the start of each year (t ≥ 1) RF(t) is calculated and the required withdrawal is made. Figure 1
tracks RF(t) over the retirement horizon for both retirees (lines are smoothed) and indicates that
Retiree 1 experienced the event Ruin(20) when funds were insufficient to make the time t=20
withdrawal. Retiree 2 did not experience a ruin event, making 30 successful withdrawals for
spending. Note that at year 10 Retiree 1 was in a better financial state than Retiree 2 (RF(10) of
0.067 vs 0.077). A Question to consider: Was an α=0.5 portfolio appropriate for Retiree 1 at all
time points during decumulation?
Figure 1
The Ruin Paths of Two Retirees
This figure details 2 ruin paths retirees may experience during decumulation and shows how the ruin
factor, RF(t), is tracked over time. We know from (3b) that, for a given α, P(Ruin(t)) is an increasing
function of RF(t-1). For a fixed t, larger RF(t) correspond to higher P(Ruin(> t)). Retiree 1 indicates it is
extremely difficult to recover once RF(t) begins a parabolic climb. RF(t) has initial value of RF(0)=WR
and we seek to keep it depressed over time. If ruin occurs RF(t) becomes negative or ∞ (not shown).

E. The Event of Financial Ruin
We define Ruin(t) as the event of financial ruin at time t. Let Ruin(≤ t) represent the
event of financial ruin occurring at or before time t and let Ruin(> t) represent the event of
financial ruin occurring after time t for t=1, 2, …, TD. (See Figure 2.)

10

Figure 2
Venn Diagram of Ruin Event
This figure expresses financial ruin as an event and demonstrates the relationship between ruin events at
different time points. For example, if Ruin(1) occurs then Ruin(≤ 2) occurs since Ruin(1) ⊆ Ruin(≤ 2).
However, Ruin(≤ 2) can occur without Ruin(1) occurring. Note that a ring of the oval represents the
event of ruin occurring at that time point. This visual helps us formulate probability statements and also
assists during the upcoming induction process when we are faced with restricted sample spaces.

It follows that,
Ruin(t)

≡

Ruin(≤ t) ∩ RuinC(≤ t-1),

(5)

which represents a single ring of the oval in Figure 2. Note that P(Ruin) can be expressed as:
1 – P(RuinC(≤ TD))

(6a)

=

1 – P(RuinC(1) ∩ RuinC(2) ∩ … ∩ RuinC(TD))

(6b)

=

1 – P(RuinC(1))*P(RuinC(2) ∩ … ∩ RuinC(TD) | RuinC(1))

(6c)

P(Ruin(> 0)) =

where the last expression holds since, for any events A, B, and C:
P(A ∩ B ∩ C) = P(A)*P(B ∩ C | A).

(7)

The retiree’s objective is to minimize P(Ruin) = P(Ruin(≤ TD)) = P(Ruin(> 0)), given WR.
As noted, the debate regarding the true value of P(Ruin) is unsettled. Further, the preferred
glide-path is subject to much disagreement. It would be helpful for the retiree to know the exact
optimal (minimum) P(Ruin) for each possible WR at every time t, along with the α required to

11

achieve this optimum. The retiree would benefit from a roadmap that informs them how their
P(Ruin) changes over time and the actions needed to maintain optimality. Ultimately such a
roadmap should be available at time t=0, guiding the retiree through decumulation.
F. Minimizing P(Ruin) at the Next Time Point
The retiree’s objective is to minimize P(Ruin) during decumulation and the choice of α at
each time point must consider the remaining horizon. A relevant sub-problem is minimizing
P(Ruin) at the next time point t ≤ TD (fixed). Assume the withdrawal at time t-1 has been made
and RF(t-1) > 0 calculated. The α that minimizes P(Ruin(t)) is found by inspecting the tail of all
ȓ(t,α) PDFs where Ruin(t) ↔ ȓ(t,α) ≤ RF(t-1). Figure 3 displays a collection of such PDFs.
Figure 3
Generic ȓ(t,α) PDFs for Various Asset Allocations (α)
This figure demonstrates how we evaluate a collection of PDFs and select one that minimizes P(Ruin(t))
recalling that Ruin(t) ↔ ȓ(t,α) ≤ RF(t-1). As we increase α the PDF of ȓ(t,α) shifts right and the variance
increases resulting in a shorter and wider distribution. Note that P(Ruin(t)) = Fȓ(t,α)(RF(t-1)), where Fx(·)
represents the CDF of the RV x. The probability of ruin is thus a left tail area under the PDF of ȓ(t,α).

In Figure 3, the density with low peak and fat tail may represent an α=1 portfolio, and the
density with high peak and smallest variance may represent the minimum volatility portfolio

12

described by Markowitz (1952). Given RuinC(t-1), the α that minimizes P(Ruin(t)) can be
determined by choosing the PDF with smallest tail area to the left of RF(t-1). This tail area is
defined by Fȓ(t,α)(RF(t-1)) where Fx(x) = P(x ≤ x) denotes the CDF of the RV x.
G. Minimizing P(Ruin) at any Time Point for Fixed TD
To minimize P(Ruin) at any time during decumulation, we use backward induction. The
induction process begins at time t=TD and each step unfolds similarly employing the expressions
in (6). Namely, we assume the retiree arrives at that time point and successfully makes their
withdrawal then recalculates RF(t) which is > 0 since RuinC(t) occurs. At t=0 the process ends
and reveals the exact optimal (minimum) P(Ruin) for all WR, along with the dynamic asset
allocation required to maintain optimality throughout decumulation. Let,
V(t, RF(t)) = Optimal (minimum) probability of ruin after time t given RF(t) > 0
α(t, RF(t)) = Value of α required to achieve V(t, RF(t)).
In terms of ruin events, the value function V(t, RF(t)) is defined using (6c) as:
V(t, RF(t)) = Min(α) P(Ruin(> t))

(8a)

= Min(α) P(Ruin(t+1) ∪ Ruin(t+2) ∪ … ∪ Ruin(TD))

(8b)

= Min(α) 1 - P(RuinC(t+1)*P(RuinC(t+2) ∩ … ∩ RuinC(TD) | RuinC(t+1)) ,

(8c)

We assume that returns from α-based portfolios at different time points are independent,
consistent with the theory of efficient markets. That is, any predictive capacity contained in past
patterns of returns is deemed already accounted for, rendering it valueless to the retiree investor.
G.1 The Induction Process for Fixed TD
Assume the retiree arrives at time t=TD and makes their final withdrawal. The restricted
sample space S={RuinC(≤ TD)} includes a single event. RF(TD) (> 0) need not be computed as
there are no future withdrawals. At t=TD, P(Ruin(>TD))=0, and V(TD, RF(TD))=0, ∀ RF(TD) > 0
serves as a boundary condition (B.C.). At any other time point t, assume the retiree makes their
withdrawal and updates RF(t) (> 0) based on the return just observed, ȓ(t,α). The retiree faces the
restricted sample space S={Ruin(t+1), Ruin(t+2), …, Ruin(TD), RuinC(≤ TD)} and seeks α to
minimize P(Ruin(> t)). Using (8c), this is equivalent to: (See Appendix C.)

13

V(t, RF(t)) =

Min

1 - (1 – Fȓ(t+1,α)(RF(t)))*
(1 – Eȓ(t+1,α)+ V t

1,

ȓ

(9)

)

,

For: 0 ≤ t ≤ TD-1, RF(t) > 0, V(TD, RF(TD)) = 0.
Optimality is achieved at α(t, RF(t)) and the expectation is wrt ȓ(t+1,α)+ = ȓ(t+1,α) | ȓ(t+1,α) > RF(t).
Assuming f(ȓ(t+1,α)) is the PDF of ȓ(t+1,α) and RF(t) is known, ȓ(t+1,α)+ has PDF:
ȓ

ȓ(t+1,α)+ = (ȓ(t+1,α) | ȓ(t+1,α) > RF(t)) ~

ȓ

,
,

ȓ

,

0

, for ȓ

,

RF t

, O.W.

(10)

1

The choice of α in (9) balances minimizing P(Ruin) at the next time point, and after the
next time point. Selecting a low volatility portfolio early in retirement to avoid ruin at the next
time point comes with a price because that portfolio results in a higher expected ruin factor,
which increases P(Ruin) after the next time point. A balance is struck at α(t, RF(t)).
H. Minimizing P(Ruin) at any Time Point for Random TD
For random TD, select SMax such that P(TD > SMax) = 0 and build discrete time hazard
probabilities P(TD=t | TD ≥ t) for t=0, 1, 2, …, SMax. Induction begins at t=SMax but now must
recognize that any time t may represent TD. The random TD value function is defined as:
VR(t, RF(t)) = Optimal (minimum) probability of ruin after time t given RF(t) > 0
αR(t, RF(t)) = Value of α required to achieve VR(t, RF(t)).
Given TD ≥ t, RuinC(t) occurs if the withdrawal at time t is successful and RuinC(≤ SMax) occurs
at time t if the withdrawal made is the last attempted. The sub-problem presented in Section II-F
must change to reflect this. We assume that TD and ȓ(t,α) are independent RVs ∀ t=1, …, SMax.
H.1 The Induction Process for Random TD
Assume the retiree arrives at time t=SMax and makes their withdrawal. RF(SMax) (> 0)
need not be calculated since there are no more withdrawals and RuinC(≤ SMax) has occurred. At
t=SMax, P(Ruin(>SMax))=0 and a B.C. is VR(SMax, RF(SMax))=0, ∀ RF(SMax) > 0. For any other
time t, assume the retiree arrives at that time point, makes their withdrawal and updates RF(t) > 0.

14

The value function VR(t, RF(t)) is defined as:
VR(t, RF(t))

= Min(α) P(Ruin(> t))

(11a)

= Min(α) 1 - P(RuinC(t+1) ∩ P(RuinC(t+2) ∩ … ∩ RuinC(SMax))

(11b)

= Min(α) P(TD > t | TD ≥ t)*[1 - {P(ȓ(t+1,α) > RF(t))*

(11c)

P(RuinC(t+2) ∩ … ∩ RuinC(SMax) | ȓ(t+1,α) > RF(t))}]
We show in Appendix D that this minimization over (0 ≤ α ≤ 1) is equivalent to:
VR(t, RF(t)) =

Min

P(TD > t | TD ≥ t)*{1 - (1 – Fȓ(t+1,α)(RF(t)))*
(1 – Eȓ(t+1,α)+ V t

1,

ȓ

,

)}

For: 0 ≤ t ≤ SMax-1, RF(t) > 0, VR(SMax, RF(SMax)) = 0.

(12)

Optimality is achieved at αR(t, RF(t)) and the expectation is over the conditional RV
ȓ(t+1,α)+.

VR(·) in (12) resembles V(·) in (9) except for the factor P(TD > t | TD ≥ t), which bounds

P(Ruin) from above when TD is random. This is intuitive since P(Ruin) at any future time point
cannot exceed the probability of living to attempt another withdrawal. We see that the fixed TD
model in (9) is a special case of (12) with P(TD=SMax)=1 or P(TD > t | TD ≥ t)=1, t=0,1,…,SMax-1
and P(TD > SMax | TD ≥ SMax)=0. V(t, RF(t)) and VR(t, RF(t)) represent value functions of dynamic
programs (DPs) with V(0, RF(0)) and VR(0, RF(0)) being the optimal (minimum) P(Ruin) at any
time point given RF(0)=WR. Time t is considered the DP’s stage, and RF(t) the state. Since the
solution finds the minimum P(Ruin) ∀ (t, RF(t)) it is also a roadmap for maintaining optimality
over time. Both V(·) and VR(·) are exact, not simulated, and optimal meaning it is not possible
to achieve a lower P(Ruin) under any other strategy using WR=RF(0) and the same stocks/bonds.
The challenge we face is solving these DPs for all RF(t) > 0, at each t. The nature of V(·)
and VR(·) will depend on the distributional assumptions made for inflation/expense-adjusted
returns. If the expression being minimized has known form then finding α(·) is handled with
calculus. Closed form expressions for V(·) and VR(·) are unlikely however. To overcome this
we discretize the α and RF(t) dimensions for 0 ≤ α ≤ 1 and RF(t) > 0, generating a numerical
approximation to the optimal solution not based on simulation.

The user can obtain an

approximation to any desired degree of accuracy by increasing the discretization precision.

15

I. Discretizing Along α and RF(t) Dimensions
Define α{·} and R{·} as the sets that result from discretization precisions Pα and PR along
the α and RF(t) dimensions, respectively with maximum RF(t) = RFMax ∀ t. That is,
α{·} = {0, 1/Pα, 2/Pα, …, 1-(1/Pα), 1}

(13a)

R{·} = {1/PR, 2/PR, 3/PR, …, RFMax-(1/PR), RFMax}.

(13b)

Along the ruin factor dimension each value (> 1/PR) in the set R{·} will serve as the midpoint of a
bucket constructed around it. Every RF(t) > 0 maps to one of the (PR)*RFMax+1 buckets, with the
last bucket containing all ruin factors > RFMax+1/(2PR). (See Figure 4.)
Figure 4
Discretization of the Positive Ruin Factor Dimension
The DPs proposed in (9) and (12) are unlikely to exhibit closed form solutions since CDFs of RVs used to
model market returns typically do not have closed functional forms. To implement these DPs we
discretize the continuous dimensions of α and RF(t). Discretizing α is trivial, but discretizing RF(t) is
more involved. In this figure we show the RF(t) discretization strategy given ruin factor precision PR and
maximum ruin factor RFMax. The discretization is contiguous (no gaps).

We will use the bucket midpoints when the discrete DP implementation requires a single
RF(t) value. Since the buckets are contiguous (no gaps), the solution provides coverage for all
RF(t) > 0. The value functions V(·) and VR(·) are derived for (t, i/PR), i=1, 2, 3, …, (PR)*RFMax,
and V(t, RF(t)) = 1 for RF(t) > RFMax + 1/(2PR). As PR→∞ the buckets collapse to the midpoint
and the discretization along the RF(t) dimension approaches the continuous solution. We handle
the conditional expectations by discretizing the probability distribution of ȓ(t+1,α)+ in a manner that
coincides with the bucket boundaries. At time t the probability that RF(t+1) falls in bucket i, i =
1, 2, …, (PR)*RFMax, (PR)*RFMax+1 given that RF(t+1) > 0 is completely determined by the CDF
of ȓ(t+1,α). The probability of falling into any bucket is precisely P(RuinC(t+1)) = P(RF(t+1) > 0)

16

which is also determined by the CDF of ȓ(t+1,α). The expectations over ȓ(t+1,α)+ in (9) and (12)
become sums over the values of V(t+1, RF(t+1)) and VR(t+1, RF(t+1)) multiplied by the
appropriate discrete conditional probability (see Appendix E). Let Vd(t, RF(t)) and αd(t, RF(t))
denote the discretized versions of V(·) and α(·) from (9), respectively. Then,
Vd(t, i/PR) =
Min

∈

1 - (1 – Fȓ(t+1,α)(i/PR))*
(1 – ∑

∗

P j ∗ V t

1, j/P

)

For 0 ≤ t ≤ TD-1, i=1, 2, 3, …, (PR)*RFMax,

(14)

w/B.C.’s Vd(TD, RF(TD)) = 0 ∀ RF(TD) > 0,
Vd(t, RF(t)) = 1 for RF(t) > RFMax+1/(2PR).
Here, P j denotes P(RF(t+1) in Bucket # j | ȓ(t+1,α) > RF(t)) with optimal ᾶ=αd(t, RF(t)).
This DP can be solved in any programming language with the user’s choice of Pα and PR. The
code operates on a grid over the (t, RF(t)) dimensions (see Figure 5). A basic (but inefficient)
strategy for coding this DP is to start at the diamond, follow the arrows, and end at the square.
Values of Vd(·) and αd(·) are derived and retained at each cell in the grid. As indicated, V(·)
from (9) exhibits structure across the t and RF(t) dimensions. Given RF(t), V(·) decreases with t,
and given t, V(·) increases with RF(t). Justifications are by indirection. Assume that one unit of
time elapses and RF(t) does not change yet V(t, RF(t)) increases. As noted, the real account
balance at time t is ($A)*RF(0)/RF(t). Therefore, the real account balance at time t+1 equals the
real account balance at time t, but P(Ruin) has increased. This is a contradiction. Similarly,
assume that at time t, RF1(t) < RF2(t) but V(t, RF1(t)) > V(t, RF2(t)). This is also a contradiction
since ($A)*RF(0)/RF1(t) > ($A)*RF(0)/RF2(t) and WR=RF(0) for both portfolios.
V(·) and VR(·) are optimal in the sense that no alternate strategy can produce a lower
P(Ruin) using the same stocks/bonds, given WR. The discrete DP is a numerical approximation
of the optimal strategy and its implementation is CPU intensive, especially for large Pα and PR.
In general, larger Pα decrease P(Ruin) while larger PR improve the approximation’s accuracy.

17

Figure 5
Discretized DP Coded as a Time by Ruin Factor Grid
The discretized version of the DP in (9) is given in (14) and it is directly extendable to the random TD
model in (12). We can solve these DPs using many strategies but generally favor solutions with shorter
processing times. One basic but inefficient strategy is shown in this figure. Here, our code would iterate
over both the time and RF(t) dimensions starting at the diamond and ending at the square. In this figure
TD is the time of final withdrawal, PR the ruin factor precision and V(·) the DP’s value function.

III. Implementations for Fixed and Random TD
We use the technique in (14) to code the DP for fixed TD in (9) with Pα=1,000,
PR=5,000, RFMax=2.75, and ER={0.5%, 0.0%}. Under this discretization 13,751 RF(t) buckets
exist at each time t. We set the unit of time to years with TD=30 resulting in a Figure 5 grid of
426,281 cells populated with (V, α). Each year the retiree calculates RF(t), slots it to a bucket
and consults the grid. The value of V(·) reflects the optimal (minimum) P(Ruin) at any future
time point up to t=TD using asset allocation α. Note that DPs require an optimal policy be
followed at the current, and all future time points reflecting the principle of optimality. To
facilitate DP coding we must make distributional assumptions for stock/bond returns.

18

A. Distributional Assumptions
We use historical annual S&P 500 total returns, 10-year Treasury Bond total returns, and
the CPI-U from 1928 – 2013 to represent stock returns, bond returns, and the inflation rate,
respectively in this analysis. We implicitly assume that future stock/bond returns originate from
the same underlying probability distributions as past returns and N=86 years of historical returns
are considered. After adjusting for inflation we find that real total stock and bond returns
originate as random samples from underlying normal distributions with (µ, σ) of (0.0825,
0.2007) and (0.0214, 0.0834), respectively. Further, a small positive correlation of ρ=0.04387
exists between these returns at each time t. (See Appendix F.)
The models in (9) and (12) make no assumption regarding the underlying distributions of
asset class returns. The only requirement is that tail probabilities for convex combinations of
correlated stock/bond returns are known or can be simulated/approximated at the bucket
boundaries over the set α{·}. The bucket boundaries are fixed and known once PR has been
selected, therefore the tail probabilities can be preprocessed. Array references would replace
CDF function calls in the code. Further, the CDF call at the bucket midpoint can be replaced by
the average boundary CDF probabilities without impacting convergence of the discrete DP to its
continuous counterpart.
In this implementation we assume stock/bond returns are iid over time, but the identical
distribution assumption is not required. Users may apply forecasts where means and variances
change as a function of time. Further, our models assume 2 asset classes but this too can be
relaxed if the CDF of linear combinations of (possibly correlated) returns across multiple asset
classes is obtainable. Clever coding may be required to keep the runtime reasonable since the
optimization would be over more than just α at each (t, RF(t)).
B. Discretized DP Results for Fixed TD
The discrete DP implementation for fixed TD described in Section III populates (V, α) for
a grid of 426,281 cells and 46 are shown in Figure 6. Here, V is the minimum P(Ruin) at any
future time point, and α is the required stock proportion. Since RF(0)=WR, the time t=0 row

19

reveals the optimal P(Ruin) for initial WR between 3.5% and 4.5%, along with the α required for
optimality during the first year. Noteworthy is that a 4% inflation-adjusted WR with ER=0.5%
yields an optimal P(Ruin) of 0.065 (success rate = 93.5%) and starting α of 41.4% in stocks.
Figure 6
Select Discrete Implementation Cells using Expense Ratio ER=0.5%
This figure displays a sample of 46 cells from the 426,281 cell optimal solution for (fixed) TD=30 years,
Pα=1,000, PR=5,000, and ER=0.5%. We start the optimal strategy at time t=0 using the asset allocation (α)
from the cell where WR=RF(0). For a retiree in standard form, the first withdrawal occurs at time t=1 and
the last at time t=TD=30. At each time t the retiree attempts their withdrawal and if successful updates
RF(t) then consults the grid for the new optimal (α). At time t, the optimal P(Ruin(> t)) is V which
changes over time. Three paths are shaded and indicate possible routes through the grid (i.e., glidepaths). When our portfolio performs well RF(t) decreases and we track leftward. When it performs
poorly RF(t) increases and we track rightward. When TD=30, the first and last α decisions are at time t=0
and t=29, respectively. The first and last withdrawals are at times t=1 and t=30, respectively.

Three paths are shaded within Figure 6 for illustration. If the retiree uses an initial 4%
WR and experiences favorable returns RF(t) decreases over time and the retiree will trace a
leftward path perhaps similar to Path 1. Stable returns close to WR result in a constant RF(t) and
the retiree’s path could resemble Path 2. Unfavorable returns increase RF(t) and the retiree could
track rightward along Path 3. As shown in Paths 1 and 2, P(Ruin) decreases in a favorable/stable
market and α declines. Thus Paths 1 and 2 roughly follow the glide-path strategy used in TD
funds that extend through the retirement date. Path 3, however reveals that the optimal equity

20

allocation roughly remains constant or increases in an unfavorable market. This is consistent
with findings by researchers who warn of poor sequences of returns early in retirement. See
Cohen et al. (2010), and Pfau and Kitces (2014). Note that a year 1 inflation/expense-adjusted
return of 1.799% would result in a rounded RF(1)=0.041 when RF(0)=0.040.
Figure 7 displays the results for ER=0.0%. The optimal P(Ruin) using a 4% initial WR is
0.043 (success rate = 95.7%) with starting α of 36.8% in equities. The same pattern of equities
decreasing during a favorable market and increasing in an unfavorable market holds true. The
95.7% success rate applies to any retiree who agrees to act optimally at all time points.
Figure 7
Select Discrete Implementation Cells using No Expense Ratio (ER=0.0%)
This figure displays the same information as Figure 6, except it is for the solution with ER=0.0%. With no
expense ratio the time t=0 ruin probabilities (V) are lower as is the asset allocation (α), which is expected.

The grids in Figures 6 and 7 apply to all retirees using any initial WR. If two retirees
begin time t=0 with different WR and end up with equal RF(t) at time t then their returns as well
as (V, α) will be identical for the remainder of retirement. This also means they will chart the
exact same course through the grid after the point of intersection, and WR plays no role in the
analysis. The intuition behind this conclusion is that the retiree with higher WR must experience
favorable returns and chart a leftward path prior to the point of intersection. While charting this

21

leftward path the retiree accumulates enough wealth to support larger withdrawals. We make no
assessment here regarding the likelihood of such intersections occurring, however.
B.1 Numerical Accuracy of the Solution
We assess the numerical approximation’s accuracy by doubling PR from 5,000 to 10,000.
This discretization creates 27,501 buckets at each time t and 852,531 total cells in the Figure 5
grid. The solution is directly comparable to Figure 7 (ER=0.0%). The larger PR changes P(Ruin)
and α modestly, but from a practical sense the results using PR=5,000 are more than adequate
when TD=30. (See Table I.)
Table I
Comparison of V and α for PR=5,000 vs. PR=10,000 with ER=0.0%
This table compares the minimum probability of ruin Vd(0, RF(0)) and the optimal asset allocation αd(0,
RF(0)) at t=0 for models with PR=5,000 vs. PR=10,000 using 3 withdrawal rates (WR). Here, PR is the
RF(t) discretization precision. Higher values of PR lead to a more accurate numerical approximation
when we discretize the DPs in (9) and (12), but also lead to longer run times. We show that very little is
gained by doubling PR from 5,000 to 10,000. The optimal solution using PR=5,000 has 426,281 total cells
and the optimal solution using PR=10,000 has 852,531 total cells when the alpha precision is Pα=1,000.

a

Time
(t)
0
0
0

RF(t)
(=WR)
0.035
0.040
0.045

Vd(0, RF(0)) [αd(0, RF(0))]
PR=5,000a
PR=10,000
0.01638 [0.314]
0.01637 [0.314]
0.04252 [0.368]
0.04251 [0.368]
0.08455 [0.445]
0.08454 [0.445]

From Figure 7, row t=0.

B.2 Consequences of Suboptimal Strategies
We assess the practical worth of the model in (9) and its implementation in Section III-B
by comparing it with suboptimal strategies. Fixed portfolios having α ∈ {0.00, 0.25, 0.50, 0.75,
1.00} are considered under the same assumptions set forth in Section III-A, with TD=30,
ER=0.0%, and WR=4%. When α is fixed, simulation can be used to estimate P(Ruin) or the DP
can be solved separately for each α with |α{·}|=1. The minimum is over one α at each cell and
V(0, 0.04) represents a numerical approximation to P(Ruin) using that strategy. The results are
evaluated against V(0, 0.04) in Figure 7. An α=0.50 portfolio generates the lowest P(Ruin)
having a success rate of 91.5%, compared with 95.7% under the optimal strategy. (See Table II.)

22

Table II
Comparable P(Ruin) using Fixed α Strategies
This table lists the P(Ruin) values using TD=30 years, ER=0.0%, and WR=4% for various fixed asset
allocation (α) strategies. These ruin probabilities are directly comparable to those in Figure 7. We can
generate these probabilities using simulation or the DP with a single α. Both methods were used and the
results are shown to be virtually identical. Our point is to show the difference in P(Ruin) between the best
fixed α strategy and the optimal strategy. The best fixed α is α=0.5 with corresponding P(Ruin) = 0.0850
(success rate = 91.5%). We saw in Figure 7 the optimal strategy yields a success rate of 95.7%. Thus,
the consequence of using a suboptimal fixed α strategy is to accept a doubling of P(Ruin).

a
b

Fixed α
0.00
0.25
0.50
0.75
1.00

Simulationa
0.3851
0.1204
0.0850
0.1070
0.1461

DP Solved for α{·}={α}b
0.3850
0.1204
0.0851
0.1070
0.1460

N=2.5 million/simulation
PR=5,000

The magnitude of difference between optimal and suboptimal P(Ruin) values increases
with TD. For example, using TD=50, ER=0.0%, and WR=4% a fixed α=0.75 generates the lowest
P(Ruin) of 0.217 (success rate = 78.3%) whereas the optimal strategy derived using PR=10,000
and Pα=1,000 generates a P(Ruin) of 0.137 (success rate = 86.3%).
B.3 Algorithm Scalability
Let k=TD, M=Pα, and N=RFMax*PR. The runtime for the discrete DP algorithm in (14)
scales at O(k*M*N2). This is because each cell in a row of the Figure 5 grid requires an expected
value computation that accesses N cells at the next time point, which occurs for M asset
allocations over all N cells in each of k rows. Therefore without pruning or parallel processing
the implementation from Section III-B.1 has a 4X longer runtime than the implementation from
Section III-B (since k and M did not change, but N is 2X larger). We can achieve significant and
scalable efficiencies in the code, however (see Appendix H).
C. Random TD Model for Multiple Retirees
Define a multi-person unit (MPU) as a group that pools retirement funds and distributes
systematic withdrawals amongst surviving members at each time point.

23

The MPU makes

withdrawals while it is alive and it is alive when at least one member is alive. Since systematic
withdrawals are used ruin can occur. The general MPU consists of K females and L males. Let
Fi and Mj represent the remaining lifetimes for female i and male j, respectively. At MPU
retirement start Fi and Mj are independent discrete RVs with know PMF and:
TD = max{F1, …, Fi, …, FK, M1, …, Mj, …, ML}

(15)

Let FTD(·) and fTD(·) represent the CDF and PMF of TD, respectively. The discrete time hazard
probabilities P(TD=t | TD ≥ t) are derived as fTD(t) / [1 - FTD(t-1)] for t=0, 1, 2, …, SMax with FTD(-1)
= 0. The PMF is constructed as fTD(t) = FTD(t) - FTD(t-1) and the CDF is built by recognizing that
the maximum of a set is less than or equal to a given value iff all members of the set are less
than or equal to that value. Namely,
TD ≤ t ↔ max{F1, …, Fi, …, FK, M1, …, Mj, …, ML} ≤ t
→

FTD(t)

(16a)

↔ (F1 ≤ t) ∩ … ∩ (FK ≤ t) ∩ (M1 ≤ t) ∩ … ∩ (MK ≤ t)

(16b)

=

(16c)

P(TD ≤ t) = P(F1 ≤ t)*…*P(FK ≤ t)*P(M1 ≤ t)*…*P(MK ≤ t).

The individual RHS probabilities in (16c) are derived from published SSA life tables and the
random TD model is solved by treating the MPU as an individual with known hazard
probabilities. The solution is an optimal decumulation strategy and the MPU can select WR by
first determining their desired success rate. Note that the same age M/F couple analysis in
Section III-C.1 is a special-case MPU having K+L=2 (K=L=1).
A risk with the random TD models we propose here is that the hazard probabilities can
become dated, rendering the solution obsolete. Suppose an MPU is formed and several members
perish unexpectedly in year 1. This MPU will have remaining hazard probabilities that differ
markedly from those used to construct the optimal strategy and the optimization should be rerun.
The same applies to a couple where one member dies early in retirement. As K+L increases the
risk lessens, meaning that couples are most exposed. They are also the easiest to adjust. Simply
transfer the living member to the random TD model solved for single retirees of that gender
starting at the appropriate time point. Random TD MPU models should thus be monitored over
time and updated when realized hazard probabilities differ from those used in the optimization.

24

C.1 Discretized DP Results for Random TD
We discretize the random TD model in (12) using the method proposed in Section III-C
for retiring same-age male/female (M/F) couples (MPU with K=L=1). Hazard probabilities are
derived from SSA life tables with time t=0 reflecting age 65. We assume the MPU is in standard
form. (Refer to Section II-B.) The first withdrawal attempt, if necessary, would occur at time
t=1 and the last at time t=SMax=48 (years). We use precisions Pα=1,000 and PR=5,000 with
ER=0.0% and WR={4%, 5%, 6%}. The results are compared against the best performing fixed α
strategy found using simulation (N=2.5 million per WR). All fixed α values between 0.0 and 1.0
were assessed in increments of 0.05. The best performing fixed α strategy using WR=4% is
α=0.45 and the corresponding P(Ruin) is 0.0421 (success rate = 95.8). The optimal strategy
starts time t=0 with α=0.356 and yields a minimum P(Ruin) of 0.0287 (success rate = 97.1%)
which is a 31.8% improvement. Comparisons for WR={4%, 5%, 6%} are shown in Table III.
Table III
P(Ruin) for Random Same-Age (65) M/F Couple TD: Optimal vs. Fixed α Strategies
In this table we present a comparison of the optimal random TD model vs. the best performing fixed α
strategy using withdrawal rates of WR={4%, 5%, 6%} for a same-age male/female couple. We assume
the retirees are in standard form with time t=0 reflecting age 65. The first withdrawal attempt occurs at
time t=1 (if TD ≥ 1) and the last at t=SMax=48 (if TD=48). The discrete time hazard probabilities P(TD=t |
TD ≥ t) were derived from life-tables at SSA.gov for t=0, 1, …, 48. The model was discretized using
precisions Pα=1,000 and PR=5,000 with ER=0.0%. The best performing fixed α solutions were found
using simulation (N=2.5 million/WR) with the test set {α: α=0.00 to 1.00 by 0.05}.

Demo

WR

M/F
Couple

4%
5%
6%

Suboptimal
P(Ruin) [Best Fixed α]
0.0421 [0.45]
0.1349 [0.60]
0.2523 [0.80]

Optimal
P(Ruin) [α at t=0]
0.0287 [0.356]
0.0978 [0.481]
0.2009 [0.672]

% Decrease in
P(Ruin)
31.8%
27.5%
20.4%

The corresponding optimal solution grid for the random TD model above is shown in
Figure 8 and reveals the familiar pattern of α decreasing when returns are favorable, and
increasing when returns are unfavorable. As noted in Section III-C, this model is exposed to
dated hazard risk and requires close monitoring with possible reoptimization over time.

25

Figure 8
Select Discrete Implementation Cells for Random TD: M/F Couple w/ER=0.0%
This figure displays a sample of 46 cells from the 673,799 cell optimal solution for the random TD model
with Pα=1,000, PR=5,000, ER=0.0%, and SMax=48. The first asset allocation decision for an MPU in
standard form is at time t=0 and the last at time t=47 (SMax-1). The first withdrawal is attempted at time
t=1 (if TD ≥ 1) and the last at time t=48 (if TD=48). Three paths are shaded as they were in Figures 6 and 7
indicating how the optimal solution changes over time under different portfolio performances. The value
of V reflects P(Ruin) at any future time point and α is the corresponding optimal asset allocation.

D. Making Adjustments over Time
We make adjustments to the optimal strategy as follows. Suppose retirees charting Path
1 in Figure 7 or 8 are dissatisfied with the decreasing α, and retirees charting Path 3 are
dissatisfied with the non-decreasing P(Ruin). Both prefer a shift to Path 2 after making the time
t=5 withdrawal. Table IV provides the real returns r(t, α) that track the RF(t) bucket midpoints for
each path. Both began time t=0 with WR=4% and an account balance of $A. In time t=0 dollars
the each withdraws (4%)*($A)/year. The Path 1 retiree can shift to Path 2 by increasing the real
withdrawal amount to (4%)*(1.143)*($A) for t >= 6. The Path 3 retiree can shift by lowering the
real withdrawal amount to (4%)*(0.889)*($A) for t >= 6.

The new rates are thus WR=

(4%)*(1.143)=4.57% and WR=(4%)*(0.889)=3.56%, respectively, if based on $A at time t=0. At
time t=6 the retirees begin their new strategy and withdraw (4.57%)*($A)* ∏i

26

1

1 I

and

(3.56%)* $A *∏i

1

1 I , respectively. Essentially, they halt their original plans and begin

anew with time t=5 as the new time t=0, and horizon of TD=25 years in the case of Figure 7. The
starting balances are (1.143)*($A)*∏i

1

1 I and (0.889)*($A)*∏i

1

1 I , respectively (see

Table IV). By shifting to Path 2, both use WR=4%=RF(5), but it is now based on their time t=5
balances. We follow the same process when provisioning for emergencies.
Table IV
Real Returns that Generate Paths 1 and 3 from Figure 7 or 8
This figure displays the rounded returns that would track the RF(t) bucket midpoints of Paths 1 and 3 in
Figures 7 and 8 (which use different α). These paths assume WR=4% and ER=0.0. At each time point we
apply the real return then subtract the real withdrawal (WR)*($A), where $A is the time t=0 account
balance. We can change WR at any time point by starting over with a new balance and consulting the grid
for the optimal α. A retiree that tracks Path 1 starts time t=5 with new balance (1.143)*($A)*∏i=1 (1+I ).
This is the new time t=0 and time t=6 reflects the new time t=1 using a modified WR.

Time
(t)
1
2
3
4
5

Real Return
r(t, α)
6.56%
6.53%
6.50%
6.48%
6.46%

Path 1
Real Account
Balance
RF(t)
.039
.038
.037
.036
(1.143)*($A) .035

Real Return
r(t, α)
1.56%
1.72%
1.87%
2.03%
2.18%

Path 3
Real Account
Balance

(0.889)*($A)

RF(t)
.041
.042
.043
.044
.045

IV. Summary and Conclusions
The discrete time models proposed in (9) and (12) yield optimal retirement decumulation
strategies. Once distributional assumptions regarding asset class returns are made these models
estimate the optimal strategy. As formulated they are intractable under common distributional
assumptions and must be discretized. The discretized solution is thus a numerical approximation
to the estimate of an optimal strategy. Since the user controls the discretization’s precision, the
approximation can be driven to any desired degree of accuracy. This leaves distributional
assumptions as the determining factor regarding how close the user’s solution is to true
optimality, and different users are sure to make different distributional assumptions.
We emphasize that these decumulation strategies are based on an objective of minimizing
the probability of ruin, not maximizing terminal wealth. Retirees seeking to maximize bequest

27

wealth may use this approach by increasing WR to gain more risk exposure, then invest unspent
funds aggressively outside of their retirement portfolio. However, it may be more prudent for
such a retiree to employ a model that was built on a wealth maximization objective, for example
that proposed by Fan, Murray, and Pittman (2013). We note that evidence suggests retirees more
often bias towards loss prevention than wealth maximization. A 2012 study by ING Retirement
Research Institute found that 80% of TD fund users prefer a portfolio that protects against
losses, and 66% of non TD fund users agreed.
A usage scenario for an advisor who seeks to employ the models proposed here is as
follows. The advisor codes or has the DP in (14) coded based on the market return assumptions
they perceive as most appropriate for the coming retirement horizon using either a fixed or
random TD, and desired levels of precision.3 The result is a grid similar in form to Figure 5. The
advisor then customizes this grid for each retiree by shading various regions the retiree indicates
they are comfortable and uncomfortable entering. A plan is then set forth at time t=0 which
indicates precisely when and what type of adjustments will be made if the retiree encroaches on a
region they have indicated is intolerable. The retiree thus takes comfort in the knowledge that a
strategy exists, and they remain fully informed of what actions will be taken, and when they will
be taken to modify that strategy based on their portfolio’s performance over time.
With this research we prove that decumulation strategies whose glide-paths are fixed at
the outset of retirement are suboptimal and this is formalized in Appendix G. Shifting from
equities to bonds reduces volatility but also the expected return and this shift is made without
regard to the retiree’s withdrawal rate. It therefore increases retirement risk when defining risk
as outliving one’s savings. When using a safe withdrawal rate strategy, the only optimal glidepath is the one that responds to market returns as they are observed over time. We have
presented a means to approximate that glide-path to any desired degree of accuracy. Lastly, we
hope this research ends the perception that the safe withdrawal rate is a simplistic rule-of-thumb
unjustified by academic theory. In fact, the theoretical framework underlying safe withdrawal
rates is rich and exhibits a natural extension to principles of optimization that are used across a
wide array of disciplines.

28

REFERENCES
Anton, Howard, 1988, Calculus with Analytic Geometry 3rd Edition (John Wiley & Sons, New
York, NY).
Bengen, William P., 1994, Determining withdrawal rates using historical data, Journal of
Financial Planning 7(4), 171-180.
Bengen, William P., 2006, Baking a withdrawal plan 'layer cake' for your retirement clients,
Journal of Financial Planning 19(8), 44-51.
Blanchett, David M., 2007, Dynamic allocation strategies for distribution portfolios:
Determining the optimal distribution glide-path, Journal of Financial Planning 20(12), 68-81.
Box, George E.P., Gwilym M. Jenkins, and Gregory C. Reinsel, 1994, Time Series Analysis
Forecasting and Control 3rd Edition (Prentice Hall, Englewood Cliffs, NJ).
Bradley, Stephen P., Arnoldo C. Hax, and Thomas L. Magnanti, 1977, Applied Mathematical
Programming (Addison-Wesley Publishing Company, Reading, MA).
Brigham, Eugene F., Michael C. Ehrhardt, 2008, Financial Management Theory and Practice
12th Edition (South-Western CENGAGE Learning, Mason, OH).
Charlson, Josh, Michael Herbst, Kailin Liu, Laura P. Lutton, and John Rekenthaler, 2009,
Target-date series research paper: Industry survey (Morningstar, Inc., Chicago, IL).
Cohen, J., Grant Gardner, and Yuan-An Fan, 2010, The date debate: Should target-date fund
glide-paths be managed "to" or "through" retirement? (Russell Research, Russell Investments,
Seattle, WA).
Cooley, Philip L., Carl M. Hubbard, and Daniel T. Walz, 1998, Retirement savings: Choosing a
withdrawal rate that is sustainable, American Association of Individual Investors Journal 10(3),
16–21.
Cooley, Philip L., Carl M. Hubbard, and Daniel T. Walz, 2011, Portfolio success rates: Where to
draw the line, Journal of Financial Planning 24(4), 48-60.
Damodaran, Aswath, Financial Database of Historical Stock and Bond Returns (New York
University, New York) <http://pages.stern.nyu.edu/~adamodar/>

29

Fan, Yuan A., Steve Murray and Sam Pittman, 2013, Optimizing Retirement Income: An
Adaptive Approach Based on Assets and Liabilities, The Journal of Retirement 1(1), 124-135.
The Federal Reserve Bank of Minneapolis, CPI-U Data 1913 – Present,
<http://www.minneapolisfed.org/community_education/teacher/calc/hist1913.cfm>
Finke, Michael, Wade D. Pfau, and David M. Blanchett, 2013, The 4 percent rule is not safe in a
low-yield world, Journal of Financial Planning 26(6), 46-55.
Guyton, Jonathan T., 2004, Decision rules and portfolio management for retirees: Is the 'safe'
initial withdrawal rate too safe? Journal of Financial Planning 17(10), 54-62.
ING Investment Management and ING Retirement Research Institute, 2012, Participant
preferences in target-date funds: Examining perceptions and expectations among target-date
users and non-users, <http://graphics8.nytimes.com/packages/pdf/tdf_white_paper.pdf>
Irlam, Gordon, and Joseph Tomlinson, 2014, Retirement income research: What can we learn
from economics? Journal of Retirement 1(4), 118-128.
Kitces, Michael E., 2012, Utility Functions in Financial Planning - A New Framework For
Decision Making? Nerd's Eye View, April 3, <http://www.kitces.com/blog/utility-functions-infinancial-planning-a-new-framework-for-decision-making/>
Kotlikoff, Laurence J., 2008, Economics' approach to financial planning. Journal of Financial
Planning 21(3), 42-52.
Markowitz, Harry, 1952, Portfolio selection, The Journal of Finance 7(1), 77–91.
Milevsky, Moshe A., and Huaxiong Huang, 2010, Spending retirement on planet vulcan: The
impact of longevity risk aversion on optimal withdrawal rates, Financial Analysts Journal 67(2),
45-58.
Milevsky, Moshe A., and Chris Robinson, 2005, A sustainable spending rate without simulation,
Financial Analysts Journal 61(6), 89-100.
Mitchell, Olivia S., and Kent Smetters, 2013, The Market for Retirement Financial Advice
(Pension Research Council) (Oxford University Press, Oxford).
Pfau, Wade D., and Michael E. Kitces, 2014, Reducing retirement risk with a rising equity glidepath. Journal of Financial Planning 27(1), 38-45.

30

Pye, Gordon B., 2000, Sustainable investment withdrawals, Journal of Portfolio Management
26(4), 73-83.
Ross, Sheldon M., 2007, Introduction to Probability Models 9th Edition (Academic Press, New
York, NY).
Scott, Jason S., William F. Sharpe, and John G. Watson, The 4% rule - At what price? Journal of
Investment Management 7(3), 31-48.
Stout, R. Gene, 2008, Stochastic optimization of retirement portfolio asset allocations and
withdrawals. Financial Services Review 17(1), 1-15.
1

These and similar figures are commonly reported in the press. See, for example, “Baby Boomers Approach 65 –
Glumly” by D. Cohn and P. Taylor of the Pew Research Center. (December 20, 2010). Web URL:
http://www.pewsocialtrends.org/2010/12/20/baby-boomers-approach-65-glumly/
2

The proof is trivial. A time t=1 withdrawal amount of (W0)*($B)*(1+I1) maintains the retiree’s purchasing power.
Replace $B with $A/(1-W0) and W0 with WR/(1+WR) then simplify. The result is (WR)*(A)*(1+I1) which is exactly
the time t=1 withdrawal amount assumed by the model’s standard form representation.

3

A full C++ implementation is provided in Appendix H.

31

October 15, 2014

©
Minimizing the Probability of Ruin in Retirement
CHRISTOPHER J. ROOK*

INTERNET APPENDIX

* The author is a consultant statistical programmer and studies in the Department of Systems Engineering
at Stevens Institute of Technology. This document accompanies the primary research paper and includes
proofs, derivations, source code, and miscellanea.

Table of Contents
Appendix A. Derivation of the Real Account Balance at Time t .................................................. 1
Appendix B. Criteria for Ruin(t) Given RuinC(t-1) ....................................................................... 2
Appendix C. Induction for Fixed TD.............................................................................................. 3
C.1 Induction at Time t=TD ....................................................................................................... 3
C.2 Induction at Time t=TD-1 .................................................................................................... 3
C.3 Induction at Time t=TD-2 .................................................................................................... 4
C.4 Induction at Time t=TD-3 .................................................................................................... 7
C.5 Induction at Time t=TD-k .................................................................................................... 9
Appendix D. Induction for Random TD ....................................................................................... 11
D.1 Induction at Time t=SMax .................................................................................................. 11
D.2 Induction at Time t=SMax-1 ............................................................................................... 11
D.3 Induction at Time t=SMax-2 ............................................................................................... 12
D.4 Induction at Time t=SMax-3 ............................................................................................... 15
D.5 Induction at Time t=SMax-k ............................................................................................... 17
Appendix E. Derivation of Ruin Factor Bucket Probabilities ..................................................... 20
Appendix F. Derivation of Historical Stock and Bond Distributions .......................................... 22
Appendix G. Miscellanea............................................................................................................. 27
G.1 The Equity Glide-Path Debate .......................................................................................... 27
G.2 Splitting the Joint PDF during Backward Induction ......................................................... 29
Appendix H. Full C++ Implementation ....................................................................................... 33

Appendix A. Derivation of the Real Account Balance at Time t
LEMMA A1: Given RuinC(≤ t), the following equality holds:
ȓ(t,α) = RF(t-1)*[1 + 1/RF(t)]

(A.1)

RF(t) = RF(t-1) / [ȓ(t,α) – RF(t-1)]

(A.2a)

→ [ȓ(t,α) – RF(t-1)]*RF(t) = RF(t-1)

(A.2b)

→ [ȓ(t,α) – RF(t-1)] = RF(t-1)/RF(t)

(A.2c)

→ ȓ(t,α) = RF(t-1)/RF(t) + RF(t-1)

(A.2d)

→ ȓ(t,α) = RF(t-1)*[1 + 1/RF(t)]

(A.2e)

Proof: By definition,

PROPOSITION A1: Given RuinC(≤ t), the real account balance at time t is ($A)*RF(0)/RF(t).
Proof: By induction we show the proposition holds for base case times t={0, 1}. Then
assuming the proposition holds at time t=t-1 we show it must also follow at time t.
Real Account Balance (t=0):

Real Account Balance (t=1):

$A

(A.3a)

= ($A)*(1)

(A.3b)

= ($A)*RF(0)/RF(0)

(A.3c)

($A)*(ȓ(1,α)) – ($A)*(WR)

(A.4a)

= ($A)*(ȓ(1,α) – WR)

(A.4b)

= ($A)*(RF(0)*[1 + 1/RF(1)] – WR)

(A.4c)

= ($A)*(RF(0)*[1 + 1/RF(1)] – RF(0))

(A.4d)

= ($A)*(RF(0)*[1 + 1/RF(1) – 1])

(A.4e)

= ($A)*RF(0)/RF(1)

(A.4f)

Real Account Balance (t=t-1): ($A)*RF(0)/RF(t-1)
Real Account Balance (t=t):

By Lemma A1.

Assume

[($A)*RF(0)/RF(t-1)]*(ȓ(t,α)) – ($A)*WR)

1

(A.5)
(A.6)

By Lemma A1.

= ($A)*[RF(0)*RF(t-1)*[1+1/RF(t)]/RF(t-1)-WR]

(A.7a)

= ($A)*[RF(0)*[1+1/RF(t)] – RF(0)]

(A.7b)

= ($A)*[RF(0)*(1+1/RF(t) – 1)]

(A.7c)

= ($A)*RF(0)/RF(t)

(A.7d)

We noted in Section II-C that the reciprocal of the ruin factor equals the number of real
withdrawals remaining. This statement is a direct result of Proposition A1. That is, 1/RF(t) = #
of real withdrawals remaining at time t, and RF(t) = 1/(# of real withdrawals remaining at time t).
Appendix B. Criteria for Ruin(t) Given RuinC(t-1)
PROPOSITION B1: Given RuinC(≤ t-1), Ruin(t) occurs iff ȓ(t,α) ≤ RF(t-1).
By Proposition A1.

Proof:
Real Account Balance (t=t-1):

($A)*RF(0)/RF(t-1)

Actual Account Balance (t=t-1):

[($A)*RF(0)/RF(t-1)]*∏i

1

1 I

(A.8b)

Actual Account Balance (t=t)1:

[($A)*RF(0)/RF(t-1)]*∏i

1

1 I

(A.8c)

(A.8a)

*(1+R(t,α))*(1 – ER)

Actual Withdrawal Amount (t=t):

($A)*(WR)*∏i

1

(A.8d)

1 I

Condition Required for Ruin(t):
≤

Account Balance (t=t)

(1+R(t,α))*(1-ER)*[($A)*RF(0)/RF(t-1)]*∏i

1

This value is pre-withdrawal.

2

1

1 I

Withdrawal Amount (t=t)

≤ ($A)*(WR)*∏i

1

1 I

(A.8e)

↔ (1+r(t,α))*(1-ER)*[1/RF(t-1)] ≤ 1

(A.8f)

↔ ȓ(t,α) ≤ RF(t-1)

(A.8g)

Appendix C. Induction for Fixed TD
C.1 Induction at Time t=TD
Assume the retiree arrives at time t=TD and makes their
last withdrawal. The restricted sample space S={RuinC(≤ TD)}
includes a single event, shown at right. The retiree need not
compute RF(TD) as there are no more withdrawals, however
RF(TD) > 0 (if computed) since RuinC(≤ TD) has occurred. At t=TD, P(Ruin(>TD))=0, and a
boundary condition (B.C.) for the value function is V(TD, RF(TD)) = 0, ∀ RF(TD) > 0.
C.2 Induction at Time t=TD – 1
Assume the retiree arrives at time t=TD-1, makes their 2nd
last withdrawal and has one remaining. The ruin factor RF(TD1) (> 0) is calculated based on the portfolio’s return just
observed, ȓ(TD-1,α). The retiree now faces the restricted sample
space S={Ruin(TD), RuinC(≤ TD)}, as shown, and seeks the α that minimizes P(Ruin(TD)). This
straight forward decision is made using the framework presented in Section II-F. Namely, the
retiree compares the tail probabilities for various asset allocations and selects the one which
minimizes P(Ruin(TD)). This optimization can be expressed as:
V(TD-1, RF(TD-1))

=

Min

P(Ruin(TD)

(C.1a)

→ V(TD-1, RF(TD-1))

=

Min

1 - P(RuinC(TD))

(C.1b)

→ V(TD-1, RF(TD-1))

=

Min

1 - P(ȓ(TD,α) > RF(TD-1))

(C.1c)

→ V(TD-1, RF(TD-1))

=

Min

1 – (1 – Fȓ(TD,α)(RF(TD-1))) ,

(C.1d)

given the known ruin factor RF(TD-1). Here, Fȓ(TD,α)(·) denotes the known/estimated CDF of
ȓ(TD,α). Note that since V(TD, RF(TD))=0, we can equivalently express V(TD-1, RF(TD-1)) as:
3

V(TD-1, RF(TD-1)) = Min

1 – (1– Fȓ(TD,α)(RF(TD-1)))*(1 – Eȓ(TD,α)+[V(TD, RF(TD))]) ,
(C.2)
Expected prob. of no ruin after
time t=TD, given RuinC(TD).

Probability of no
ruin at time t=TD.

with optimal ᾶ=α(TD-1, RF(TD-1)).2
C.3 Induction at Time t=TD – 2
Assume the retiree arrives at time t=TD-2, makes their
3rd last withdrawal and has two remaining. The ruin factor
RF(TD-2) (> 0) is calculated based on the portfolio’s return just
observed, ȓ(TD-2,α). The retiree now faces the restricted sample
space S={Ruin(TD-1), Ruin(TD), RuinC(≤ TD)}, as shown, and seeks to make the optimal asset
allocation decision to minimize P(Ruin(>TD-2)) = P(Ruin(TD-1) ∪ Ruin(TD)), which is the
probability of ruin at any future time point. Using (6b) and (6c), we express P(Ruin(>TD-2)) as:
P(Ruin(TD-1) ∪ Ruin(TD))

= 1 - P(RuinC(TD-1) ∩ RuinC(TD))

(C.3a)

= 1 - P(RuinC(TD-1))*P(RuinC(TD) | RuinC(TD-1)).

(C.3b)

The optimal value for this probability was just derived at time t=TD-1 for all RF(t) > 0.
The value function is expressed as:
V(TD-2, RF(TD-2)) = Min
→ V(TD-2, RF(TD-2)) = Min

1 - P(RuinC(TD-1))*P(RuinC(TD) | RuinC(TD-1)) (C.4a)
1 - P(ȓ(TD-1,α) > RF(TD-2))*
(C.4b)
P(ȓ(TD,ᾶ) > RF(TD-1) | ȓ(TD-1,α) > RF(TD-2))

→ V(TD-2, RF(TD-2))3 = Min

1 - P(ȓ(TD-1,α) > RF(TD-2))*
ȓ

∩ ȓ

,ᾶ

ȓ

2

,

(C.4c)

,

The right-most term is added for notational convenience. Recall that V(TD, RF(TD)) = 0 ∀ RF(TD) > 0 and it is
therefore the expected value of zero. Further, ȓ(TD,α)+ = (ȓ(TD,α) | ȓ(TD,α) > RF(TD-1)) as given in Section II-G.1.
3
Our convention is to let ᾶ refer to an alpha that is optimal at a future time point, and let α represent one that is
being optimized over at the current time point. The optimal ᾶ is always required to minimize the current V(·).

4

The numerator of the ratio in (C.4c) reflects the probability that ruin is avoided at both
times t=TD-1 and TD. By definition, this is the volume of the joint PDF f(ȓ(TD-1,α), ȓ(TD,ᾶ)) over the
region defined in the probability statement. This region exists in the ȓ(TD-1,α) – ȓ(TD,ᾶ) plane and the
joint PDF defines a 3-dimensional object that rests on the plane. We calculate the required
volume by integrating the joint PDF over the given region. Since the integration limits for ȓ(TD,ᾶ)
depend on ȓ(TD-1,α) we must handle ȓ(TD,ᾶ) first, where ȓ(TD,ᾶ) ranges from RF(TD-1) to ∞. At this
induction step ȓ(TD-1,α) ranges from the constant RF(TD-2) to ∞. This region of the ȓ(TD-1,α) – ȓ(TD,ᾶ)
plane is shown below as the cross-hatched section of Figure A1. Assuming
returns are bell-shaped, then a 3-dimensional hill object (see right) depicts
f(ȓ(TD-1,α), ȓ(TD,ᾶ)). The required probability is the volume of this object over
the region shown, and we must evaluate it across all α at time t=TD-2. Varying α changes the
shape and location of the hill, and thus the probability. The denominator of the ratio in (C.4c) is
the volume of the same solid to the right of RF(TD-2). We seek α that minimizes the entire
expression in (C.4c) where these probabilities are two components, see (C.4d) below.
Figure A1. Region under f(ȓ(TD-1,α), ȓ(TD,ᾶ)) for Derivations
of the Numerator and Denominator in (C.4c)

5

→ V(TD-2, RF(TD-2)) =
Min

Since RF(TD-1) is a function of ȓ(TD-1,α), these integrals must
remain nested and the ordering cannot be interchanged.
1 - (1 – Fȓ(TD-1,α)(RF(TD-2)))*
ȓ

,ᾶ

,ȓ

ȓ

,

ȓ

,ᾶ

ȓ

,

ȓ

(C.4d)

,

,

By conditioning and assuming returns are independent between time points we split the joint
PDF in the numerator of the ratio (see Appendix G.2):

→ V(TD-2, RF(TD-2)) =
Min

This term is [1 – V(TD-1, RF(TD-1))] found earlier. An
optimal policy must be followed at each stage of the DP.
1 - (1 – Fȓ(TD-1,α)(RF(TD-2)))*
ȓ

ȓ

,

ȓ

ȓ

,

,ᾶ

ȓ

,

(C.4e)

,

By definition, this integral is the expected
value of [1 – V(·)] over the R.V. ȓ(TD-1,α)+.

→ V(TD-2, RF(TD-2)) =
Min

ȓ

,ᾶ

1 - (1 – Fȓ(TD-1,α)(RF(TD-2)))*
ȓ

,

,

ȓ

,

ȓ
ȓ

,

(C.4f)

,

Now, since RF(TD-1) is a function of ȓ(TD-1,α), namely,
RF(TD-1)

=

ȓ

,

(C.5)

the value function in (C.4f) can be written as4,

4

An intuitive explanation of the need to take expectations is that V(TD-1, X), was already found for all positive ruin
factors X, treating X as constant. At the current induction step, it is discovered that X is random with known PDF
under our control via α. In the optimization over α, EX[V(TD-1, X)] is then evaluated across the various PDFs of X.

6

V(TD-2, RF(TD-2)) =
Min

RF(TD-1)
1 - (1 – Fȓ(TD-1,α)(RF(TD-2)))*
(1 – Eȓ(TD-1,α)+ V T

1,

ȓ

)

,

(C.6)

with optimality achieved at ᾶ=α(TD-2, RF(TD-2)) and the expectation over the conditional RV
ȓ(TD-1,α)+ = (ȓ(TD-1,α) | ȓ(TD-1,α) > RF(TD-2)) where {ȓ(TD-1,α) > RF(TD-2)} ≡ {RF(TD-1) > 0}.
C.4 Induction at Time t=TD – 3
Assume the retiree arrives at time t=TD-3, makes their
4th-last withdrawal and has 3 remaining. 5 The ruin factor
RF(TD-3) (> 0) is calculated based on the portfolio return just
observed, ȓ(TD-3,α). The retiree is facing the restricted sample
space S={Ruin(TD-2), Ruin(TD-1), Ruin(TD), RuinC(≤ TD)} (shown at right) and seeks to make
the optimal asset allocation decision to minimize P(Ruin(TD-2) ∪ Ruin(TD-1) ∪ Ruin(TD)),
which is the probability of ruin at any future time point. Under the restricted sample space,
P(Ruin(>TD-3)) is now defined as:
P(Ruin(TD-2) ∪ Ruin(TD-1) ∪ Ruin(TD))
= 1 - P(RuinC(TD-2) ∩ RuinC(TD-1) ∩ RuinC(TD))

(C.7a)

= 1 - P(RuinC(TD-2))*P(RuinC(TD-1) ∩ RuinC(TD) | RuinC(TD-2))

(C.7b)

Note that the optimal value for this probability was derived at time t=TD-2 for all RF(t) > 0.
The value function is thus:
V(TD-3, RF(TD-3)) = Min

1 - P(RuinC(TD-2))*
C

C

C

(C.8a)

P(Ruin (TD-1) ∩ Ruin (TD) | Ruin (TD-2))
5

Induction at time t=TD-3 is nearly identical to induction at time t=TD-2, and the process is generalized for time
t=TD-k next, then reported for any time t in Section II-G.1.

7

→ V(TD-3, RF(TD-3)) =
Min

1 - P(ȓ(TD-2,α) > RF(TD-3))*
(C.8b)

P(ȓ(TD-1,ᾶ) > RF(TD-2) ∩ ȓ(TD,ᾶ) > RF(TD-1) | ȓ(TD-2,α) > RF(TD-3))

We require an optimal policy be followed at each future stage, and these ᾶ
reflect those optimal values. This V(·) cannot take a minimum value otherwise.
Applying the definition of
a conditional probability.

→ V(TD-3, RF(TD-3)) =
Min

1 - P(ȓ(TD-2,α) > RF(TD-3))*
ȓ

,ᾶ

∩ȓ

,ᾶ

ȓ

(C.8d)

1 - (1 – Fȓ(TD-2,α)(RF(TD-3)))*
ȓ

,ᾶ

,ȓ

,ᾶ

ȓ

,ȓ

ȓ

,

ȓ

,

ȓ

,ᾶ

ȓ

,ᾶ

,

,

By conditioning and independence we
split the joint PDF (see Appendix G.2).

→ V(TD-3, RF(TD-3)) =
Min

(C.8c)

,

The multivariate density of the next 3 real
returns integrated over the condition of no ruin.

→ V(TD-3, RF(TD-3)) =
Min

∩ ȓ

,

(C.8e)

1 - (1 – Fȓ(TD-2,α)(RF(TD-3)))*
ȓ

ȓ

,

ȓ

→ V(TD-3, RF(TD-3)) =
Min

,ᾶ

,ȓ

ȓ

,ᾶ

ȓ

,

ȓ

,ᾶ

ȓ

,ᾶ

,

Note that the retiree has some control over the
next ruin factor via choice of α at this time point.
1 - (1 – Fȓ(TD-2,α)(RF(TD-3)))*
ȓ

,

,

ȓ

8

,

,

ȓ
ȓ

,

,

(C.8f)

Since RF(TD-2) is a function of ȓ(TD-2,α), namely,
RF(TD-2)

=

ȓ

(C.8g)
,

the expression above is by definition the expected value of [1 - V(TD-2, RF(TD-2))] over the
conditional RV ȓ(TD-2,α)+, and the value function can be written as:
V(TD-3, RF(TD-3)) =
Min

RF(TD-2)
1 - (1 – Fȓ(TD-2,α)(RF(TD-3)))*
(C.9)

+

(1 – Eȓ(TD-2,α) V T

2,

ȓ

,

)

Optimality is achieved at ᾶ=α(TD-3, RF(TD-3)) and the expectation is over the conditional RV
ȓ(TD-2,α)+ = (ȓ(TD-2,α) | ȓ(TD-2,α) > RF(TD-3)) where {ȓ(TD-2,α) > RF(TD-3)} ≡ {RF(TD-2) > 0}.
C.5 Induction at Time t=TD – k
Assume the retiree arrives at time t=TD-k, for k=0,1, …,
TD-1 and makes their (k+1)th-last withdrawal with k remaining.
The ruin factor RF(TD-k) (> 0) is calculated based on the
portfolio return just observed, ȓ(TD-k,α). The retiree is facing the
restricted sample space S={Ruin(TD-k+1), Ruin(TD-k+2), …, Ruin(TD), RuinC(≤ TD)} (shown at
right) and seeks the optimal asset allocation to minimize P(Ruin(>TD-k)) = P(Ruin(TD-k+1) ∪
Ruin(TD-k+2) ∪ … ∪ Ruin(TD)), which is the probability of ruin at any future time point. Under
the restricted sample space we express P(Ruin(>TD-k)) as:
P(Ruin(TD-k+1) ∪ Ruin(TD-k+2) ∪ … ∪ Ruin(TD))
= 1 - P(RuinC(TD-k+1) ∩ RuinC(TD-k+2 ) ∩ … ∩ RuinC(TD))

(C.10a)

= 1 - P(RuinC(TD-k+1))*P(RuinC(TD-k+2) ∩ … ∩ RuinC(TD) | RuinC(TD-k+1))

(C.10b)

The optimal value for this probability is derived at time t=TD-k+1 for all RF(t) > 0.
9

The value function is given by:
V(TD-k, RF(TD-k)) =
1 - P(RuinC(TD-k+1))*

Min

C

C

(C.10c)

C

P(Ruin (TD-k+2) ∩ … ∩ Ruin (TD) | Ruin (TD-k+1))
Let,

ȓ
Ȓ

,

=

,ᾶ

ȓ
ȓ

⋮

,ᾶ

, and,

:⋂

,

ȓ

RF t

,ᾶ

1

(C.10d)

.

,ᾶ

At time t=TD-k the vector Ȓ(TD-k+2,TD) will hold the random returns at all time points after t=TDk+1, assuming optimal asset allocations are used. The set

will represent the space

,

in k-2 dimensions over which Ȓ(TD-k+2,TD) satisfies the condition of RuinC(>TD-k+1).
→ V(TD-k, RF(TD-k)) =
Min

1 - P(ȓ(TD-k+1, α) > RF(TD-k))*
P(Ȓ

,

∈

,

1 - P(ȓ(TD-k+1, α) > RF(TD-k))*
ȓ

,

∩ Ȓ

→ V(TD-k, RF(TD-k)) =

∈

,

(C.10f)

,

,

ȓ

Min

(C.10e)

Applying the definition of
a conditional probability.

→ V(TD-k, RF(TD-k)) =
Min

| ȓ(TD-k+1, α) > RF(TD-k))

The multivariate density of the next k real
returns integrated over the condition of no ruin.

(C.10g)

1 - (1 – Fȓ(TD-k+1,α)(RF(TD-k)))*
,

ȓ

,

ȓ

10

,Ȓ

Ȓ

,
,

ȓ

,
,

ȓ

,

→ V(TD-k, RF(TD-k)) =

(C.10h)

By conditioning and independence
(see Appendix G.2).

Min

1 - (1 – Fȓ(TD-k+1,α)(RF(TD-k)))*
ȓ

,

Ȓ

,

ȓ

ȓ

,

ȓ

,

,

,

By definition, the bracketed term
in (C.10h) is P(RuinC(>TD-k+1)).

→ V(TD-k, RF(TD-k)) =
Min

Ȓ

,

(C.10i)

1 - (1 – Fȓ(TD-k+1,α)(RF(TD-k)))*
ȓ

,

,

ȓ

ȓ
ȓ

,

,

,

→ V(TD-k, RF(TD-k)) =
Min

RF(TD-k+1)

1 - (1 – Fȓ(TD-k+1,α)(RF(TD-k)))*

(C.10j)
+

(1 – Eȓ(TD-k+1,α) V T

k

1,

ȓ

,

)

Optimality is achieved at ᾶ=α(TD-k, RF(TD-k)) and the expectation is over the conditional RV
ȓ(TD-k+1,α)+ = (ȓ(TD-k+1,α) | ȓ(TD-k+1,α) > RF(TD-k)) where {ȓ(TD-k+1,α) > RF(TD-k)} ≡ {RF(TD-k+1) > 0}.
Appendix D. Induction for Random TD
D.1 Induction at Time t=SMax
Assume the retiree arrives at time t=SMax and makes their last withdrawal. RF(SMax) (> 0)
need not be calculated since there are no more withdrawals and RuinC(≤ SMax) has occurred. At
time t=SMax, P(Ruin(>SMax))=0 and a B.C. for the value function is VR(SMax, RF(SMax))=0, ∀
RF(SMax) > 0.
D.2 Induction at Time t=SMax – 1
Assume the retiree arrives at time t=SMax-1, makes their withdrawal and has at most 1
remaining. RF(SMax-1) (> 0) is calculated based on the portfolio’s return just observed, ȓ(SMax-1,α).
11

The retiree seeks α to minimize P(Ruin(>SMax-1)) = P(Ruin(SMax)), which we express as:
VR(SMax-1, RF(SMax-1)) = Min
→ VR(SMax-1, RF(SMax-1)) = Min

(D.1a)

1 – P(RuinC(SMax))

(D.1b)

RuinC(SMax) ≡ (death before time t=SMax withdrawal attempt)
∪ (live until t=SMax and successfully make withdrawal)

→ VR(SMax-1, RF(SMax-1)) =
Min

P(Ruin(SMax))

1 – [P(TD=SMax-1 | TD ≥ SMax-1)*(1)

(D.1c)

+ P(TD > SMax-1 | TD ≥ SMax-1)*P(ȓ(SMax,α) > RF(SMax-1))]
This term =
P(TD = SMax | TD ≥ SMax-1)
→ VR(SMax-1, RF(SMax-1)) =
Min

P(TD > SMax-1 | TD ≥ SMax-1)*[1 - P(ȓ(SMax,α) > RF(SMax-1))]

(D.1d)

→ VR(SMax-1, RF(SMax-1)) =
Min

P(TD > SMax-1 | TD ≥ SMax-1)*[1 – (1 – Fȓ(SMax,α)(RF(SMax-1)))] , (D.1e)

given the known ruin factor RF(SMax-1). We can alternatively express VR(SMax-1, RF(SMax-1)) as:
VR(SMax-1, RF(SMax-1)) =
Min

P(TD > SMax-1 | TD ≥ SMax-1)*

RF(SMax) =
RF(SMax-1)/(ȓ(SMax,α) - RF(SMax-1))
(D.1f)

[1 – (1 – Fȓ(SMax,α)(RF(SMax-1)))*(1 – Eȓ(SMax,α)+[VR(SMax, RF(SMax))])]
P(RuinC(SMax)) given
TD > SMax -1.

Expected prob. of no ruin after time t=SMax,
given RuinC(SMax). (This expression = 1.)

with optimal ᾶ=αR(SMax-1, RF(SMax-1)).
D.3 Induction at Time t=SMax – 2
Assume the retiree arrives at time t=SMax-2, makes their withdrawal and has at most 2
remaining. RF(SMax-2) (> 0) is calculated based on the portfolio return just observed, ȓ(SMax-2,α).
12

The retiree seeks α to minimize P(Ruin(>SMax-2)) = P(Ruin(SMax-1) ∪ Ruin(SMax)), which for a
given α can be expressed as:
P(Ruin(>SMax-2))
= P(Ruin(SMax-1) ∪ Ruin(SMax))

(D.2a)

= 1 - P(RuinC(SMax-1) ∩ RuinC(SMax))

(D.2b)

= 1 – P(TD=SMax-2 | TD ≥ SMax-2)*(1)
+ P(TD > SMax-2 | TD ≥ SMax-2)*P(ȓ(SMax-1,α) > RF(SMax-2) ∩ RuinC(SMax))

(D.2c)

= P(TD > SMax-2 | TD ≥ SMax-2)*
[1 - P(ȓ(SMax-1,α) > RF(SMax-2))*P(RuinC(SMax) | ȓ(SMax-1,α) > RF(SMax-2)) ]

(D.2d)

P(RuinC(SMax) | ȓ(SMax-1,α) > RF(SMax-2))
=

ȓ

,

*P(ȓ(SMax-1,α)

> RF(SMax-2) ∩ RuinC(SMax))

P(ȓ(SMax-1,α) > RF(SMax-2) ∩ RuinC(SMax))

(D.2e)

These 2 events are independent since
ȓ(SMax-1,α) and TD are independent R.V.s.

= P[ (ȓ(SMax-1,α) > RF(SMax-2) ∩ (TD=SMax-1 | TD > SMax-2))
∪

These 2 events are
mutually exclusive.

(ȓ(SMax-1,α) > RF(SMax-2) ∩ (TD=SMax | TD > SMax-2) ∩ ȓ(SMax,ᾶ) > RF(SMax-1)) ]

(D.2f)

These 2 events are NOT independent since
RF(SMax-1) is a function of ȓ(SMax-1,α).
The event of making the withdrawal at SMax-1 and also avoiding ruin at SMax can happen 2 ways:
Make the withdrawal at SMax-1 and experience death before SMax, or make the withdrawal at
SMax-1, then again at SMax. Note that all this unfolds given TD > SMax-2 ≡ TD ≥ SMax-1, see (D.2c).
The intersection of events over TD and (ȓ(SMax-1,α), ȓ(SMax,ᾶ)) are independent since random
market returns are independent of the retiree’s random time of final withdrawal. Also, the union
of events requiring TD=SMax-1 and TD=SMax are mutually exclusive since both cannot occur

13

simultaneously. Finally, the events ȓ(SMax-1,α) > RF(SMax-2) and ȓ(SMax,ᾶ) > RF(SMax-1) are not
independent since RF(SMax-1) = RF(SMax-2) / [ȓ(SMax-1,α) - RF(SMax-2)]. We thus express (D.2f) as:
P(ȓ(SMax-1,α) > RF(SMax-2) ∩ RuinC(SMax))
= P(ȓ(SMax-1,α) > RF(SMax-2))*P(TD=SMax-1 | TD ≥ SMax-1)

(D.2g)

+ P(ȓ(SMax-1,α) > RF(SMax-2) ∩ ȓ(SMax,ᾶ) > RF(SMax-1))*P(TD=SMax | TD ≥ SMax-1)
= P(TD=SMax-1 | TD ≥ SMax-1)*

(D.2h)

ȓ

ȓ

,

,

The joint PDF is split (see Appendix G.2).

+ P(TD=SMax | TD ≥ SMax-1)*
ȓ
=

ȓ

,

,

,ȓ

ȓ

,ᾶ

∗

ȓ

,ᾶ

ȓ

,

(D.2i)

,

ȓ

P(TD=SMax-1 | TD ≥ SMax-1) + P(TD=SMax | TD ≥ SMax-1)*

,ᾶ

ȓ

,ᾶ

If the optimal α=ᾶ is used at t=SMax-1 then this term is precisely Maxα[P(RuinC(SMax))]
given the retiree arrives at t=SMax-1 and successfully makes the withdrawal, which by
definition is 1 – VR(SMax-1, RF(SMax-1)), see (D.1c). It will be shown below that α=ᾶ
must be used at t=SMax-1 to minimize VR(SMax-2, RF(SMax-2)) at t=SMax-2.
Substituting the term from (D.2i) back into (D.2e) reveals it is nothing more than the expectation
of [1 – VR(SMax-1, RF(SMax-1))] over the conditional RV ȓ(SMax-1,α)+. If using α=ᾶ at t=SMax-1, we
can express the original probability from (D.2a) for a given α at t=SMax-2 with RF(SMax-2) as:
P(Ruin(>SMax-2))
RF(SMax-1)

= P(TD > SMax-2 | TD ≥ SMax-2)*
[1- P(ȓ(SMax-1,α) > RF(SMax-2))*(1- Eȓ(SMax-1,α)+ V S

1,

ȓ

,

)] (D.2j)

For P(Ruin(>SMax-2)) to be minimized over all α at t=SMax-2, the function in [·] must take
its minimum value at each RF(SMax-1), which occurs precisely at VR(SMax-1, RF(SMax-1)).
14

Finally, we express the value function at t=SMax-2 as:
VR(SMax-2, RF(SMax-2)) =
Min

(D.2k)

P(TD > SMax-2 | TD ≥ SMax-2)*
{1- (1 - Fȓ(SMax-1,α)(RF(SMax-2)))*(1 - Eȓ(SMax-1,α)+[VR(SMax-1, RF(SMax-1))])}

with optimality achieved at ᾶ=αR(SMax-2, RF(SMax-2)).
D.4 Induction at Time t=SMax – 3
Assume the retiree arrives at time t=SMax-3, makes their withdrawal (at most 3 remain),
and updates RF(SMax-3) (> 0) based on the portfolio return just observed, ȓ(SMax-3,α). The retiree
seeks α to minimize P(Ruin(>SMax-3)) = P(Ruin(SMax-2) ∪ P(Ruin(SMax-1) ∪ Ruin(SMax)), which
for a given α can be expressed as:
P(Ruin(>SMax-3))
= P(Ruin(SMax-2) ∪ Ruin(SMax-1) ∪ Ruin(SMax))

(D.3a)

= 1 - P(RuinC(SMax-2) ∩ RuinC(SMax-1) ∩ RuinC(SMax))

(D.3b)

=1–

P(TD=SMax-3 | TD ≥ SMax-3)*(1)
+ P(TD > SMax-3 | TD ≥ SMax-3)*

(D.3c)

P(ȓ(SMax-2,α) > RF(SMax-3) ∩ RuinC(SMax-1) ∩ RuinC(SMax))
(D.3d)

= P(TD > SMax-3 | TD ≥ SMax-3)*

[1 - P(ȓ(SMax-2,α) > RF(SMax-3))*P(RuinC(SMax-1) ∩ RuinC(SMax) | ȓ(SMax-2,α) > RF(SMax-3)) ]

P(RuinC(SMax-1) ∩ RuinC(SMax) | ȓ(SMax-2,α) > RF(SMax-3))
=

ȓ

,

*P(ȓ(SMax-2,α)

> RF(SMax-3) ∩ RuinC(SMax-1) ∩ RuinC(SMax))

15

(D.3e)

Given TD > SMax-3, there
are 3 mutually exclusive
ways this intersection of
3 events can occur:

P(ȓ(SMax-2,α) > RF(SMax-3) ∩ RuinC(SMax-1) ∩ RuinC(SMax))
= P[((ȓ(SMax-2,α)>RF(SMax-3) ∩ (TD=SMax-2 | TD ≥ SMax-2))
∪

(ȓ(SMax-2,α)>RF(SMax-3) ∩ (TD=SMax-1 | TD ≥ SMax-2) ∩ ȓ(SMax-1,ᾶ)>RF(SMax-2))

(D.3f)

∪
(ȓ(SMax-2,α) >RF(SMax-3) ∩ (TD=SMax | TD ≥ SMax-2) ∩ ȓ(SMax-1,ᾶ)>RF(SMax-2) ∩ ȓ(SMax,ᾶ)>RF(SMax-1))
Invoking the independence of TD and ȓ(·), (D.3f) can be written as:
= P(TD=SMax-2 | TD ≥ SMax-2)*P(ȓ(SMax-2,α)>RF(SMax-3)
(D.3g)

+ P(TD=SMax-1 | TD ≥ SMax-2)*P(ȓ(SMax-2,α)>RF(SMax-3) ∩ ȓ(SMax-1,ᾶ)>RF(SMax-2))

+ P(TD=SMax | TD ≥ SMax-2)*P(ȓ(SMax-2,α)>RF(SMax-3) ∩ ȓ(SMax-1,ᾶ)>RF(SMax-2) ∩ ȓ(SMax,ᾶ)>RF(SMax-1))
= P(TD=SMax-2 | TD ≥ SMax-2)*
ȓ

ȓ

,

,

We split the joint PDF (see Appendix G.2).

+ P(TD=SMax-1 | TD ≥ SMax-2)*
ȓ

,

,ȓ

ȓ

,ᾶ

ȓ

,ᾶ

(D.3h)

,

We split the joint PDF (see Appendix G.2).

+ P(TD=SMax | TD ≥ SMax-2)*
ȓ

,ᾶ

,ȓ

,ᾶ

,ȓ

ȓ

,

,ᾶ

ȓ

,ᾶ

ȓ

,

Since the sum of integrals with respect to ȓ(SMax-2,α) is the integral of the sum, (D.3h) becomes:
=

ȓ

,

∗

ȓ

(D.3i)

,

P(TD = SMax-2 | TD ≥ SMax-2)
ȓ

+ P(TD = SMax-1 | TD ≥ SMax-2)*

ȓ

,ᾶ

,ᾶ

+ P(TD = SMax | TD ≥ SMax-2)*
ȓ

,ᾶ

,ȓ

16

,ᾶ

ȓ

,ᾶ

ȓ

,ᾶ

This entire term is recognized as Maxα[1 - P(Ruin(>SMax-2))], (see D.2c) given that
the retiree arrives at t=SMax-2 and successfully makes the withdrawal, which is by
definition 1–VR(SMax-2, RF(SMax-2)). Following an optimal policy at future time
points is required to achieve a minimum for VR(SMax-3, RF(SMax-3)) at t= SMax-3.
Substituting (D.3i) back into the RHS of (D.3e) reveals it is the expected value of 1–VR(SMax-2,
RF(SMax-2)) with respect to the RV ȓ(SMax-2,α)+. This leaves (D.3d), assuming we use optimal
decumulation policies at all future time points, as:
P(Ruin(>SMax-3))
RF(SMax-2)

= P(TD > SMax-3 | TD ≥ SMax-3)*
[1- P(ȓ(SMax-2,α) > RF(SMax-3))*(1- Eȓ(SMax-2,α)+ V S

2,

ȓ

,

)] (D.3j)

The value function VR(SMax-3, RF(SMax-3)) is the minimum of (D.3j) over 0 ≤ α ≤ 1, namely:
VR(SMax-3, RF(SMax-3)) =
Min

(D.3k)

P(TD > SMax-3 | TD ≥ SMax-3)*
{1 - (1 - Fȓ(SMax-2,α)(RF(SMax-3)))*(1 - Eȓ(SMax-2,α)+[VR(SMax-2, RF(SMax-2))])}

with optimality achieved at ᾶ=αR(SMax-3, RF(SMax-3)).
D.5 Induction at Time t=SMax – k
Assume the retiree arrives at time t=SMax-k, makes their withdrawal (at most k remain),
and updates RF(SMax-k) (> 0) based on the portfolio return just observed, ȓ(SMax-k,α). The retiree
seeks α to minimize P(Ruin(>SMax-k))=P(Ruin(SMax-k+1) ∪ P(Ruin(SMax-k+2) ∪ … ∪ Ruin(SMax)),
which for a given α=αR(SMax-k, RF(SMax-k)) can be expressed as:
P(Ruin(>SMax-k))
=P ⋃
= 1 - P( ⋂

(D.4a)

Ruin t
Ruin t )

(D.4b)
17

=1–

P(TD=SMax-k | TD ≥ SMax-k)*(1)
+ P(TD > SMax-k | TD ≥ SMax-k)*

(D.4c)
Ruin t )

P(ȓ(SMax-k+1,α) > RF(SMax-k) ∩ [ ⋂

= P(TD > SMax-k | TD ≥ SMax-k)*{1 - P(ȓ(SMax-k+1,α) > RF(SMax-k))*
Ruin t | ȓ(SMax-k+1,α) > RF(SMax-k)) }

P(⋂

ȓ

⋃

k 1,

Ruin t

∗ P( ⋂

,

P{ ȓ SMax

(D.4d)

RF SMax
T

k ∩

T

t |T

S

S

k
k

(D.4e)

∩ ȓ(SMax-k+1,α) > RF(SMax-k))

1|T

S

1 ∩ ⋂

k
ȓ

This is: Ȓ

1 ∪
RF t

t2 ,ᾶ

,

(D.4f)

∈

1

}
,

.

In (D.4e) we are operating under the assumption that (TD > SMax-k | TD ≥ SMax-k), see (D.4c).
Having made it to time t=SMax-k+1 we avoid ruin at time t=SMax-k+1 iff [ȓ(SMax-k+1,α) > RF(SMax-k)],
and we avoid ruin at all times t > SMax-k+1 iff ⋂

Ruin t , which is equivalent to:

(TD = SMax-k+1 ∪ [(TD > SMax-k+1 | TD ≥ SMax-k+1) ∩ (RuinC(>SMax-k+1))]).

(D.4g)

This explains the decomposition of (D.4e) to (D.4f), where (D.4f) also expands the bracketed [·]
event to all possible remaining values of TD. These separate TD events form the basis of the
union operator “∪” and are mutually exclusive thus can be summed. The events over the
intersection operator “∩” cannot be multiplied because they are not independent (e.g., see
Figure A1). For any events A, B, and C, the distributive property of the set operator “∩” is such
that A ∩ (B ∪ C) ≡ (A ∩ B) ∪ (A ∩ C). Using this property we distribute the event [ȓ(SMax-k+1,α) >
RF(SMax-k)] over all events separated by “∪” in (D.4f) allowing it to be written as:

18

T

P{

S

1|T

k

T

⋃

S

k

t |T

1 ∩ ȓ

,α

k

1 ∩ ⋂

S

RF S

k
ȓ

t2 ,ᾶ

(D.4h)

∪
RF t

1

}

We have seen this quantity in previous induction steps. Borrowing the notation from (C.10d)
that introduced Ȓ

and

,

, and using the results from Appendix G.2 to move the

,

integration over RV ȓ(SMax-k+1,α) outside of the sum, we decompose (D.4h) as:
ȓ
∑

,

PT

* PT
t |T

S

S

1|T

k
k

1 ∗

S

,

k
Ȓ

(D.4i)

1 +
,

Ȓ

,

ȓ

,

The quantity in (D.4i) is divided by P(ȓ(SMax-k+1,α) > RF(SMax-k)) in (D.4e) leaving it as the
expected value of [·] over the RV ȓ(SMax-k+1,α)+. Assuming an optimal policy is followed at all t >
SMax-k+1 we recognize [·] as the maximum P(RuinC(>SMax-k+1)) given the retiree arrives at time
t=SMax-k+1 and makes their withdrawal, which is precisely 1 – VR[SMax-k+1, RF(SMax-k+1)].
Since our objective at t=SMax-k is to minimize P(Ruin(>SMax-k)) we must use the maximum
value of [·] ∀ RF(SMax-k). (This follows since for any 2 functions h(x)>0 and g(x)>0, h(x) ≥ g(x)
∀ x → EX[h(x)] ≥ EX[g(x)]. The proof can proceed by contradiction and is trivial since density
functions cannot take negative values.) Finally, we minimize (D.4d) over (0 ≤ α ≤1) expressing
the value function as:
VR(SMax-k, RF(SMax-k)) =
Min

(D.4j)

P(TD > SMax-k | TD ≥ SMax-k)*{1 - (1 - Fȓ(SMax-k+1,α)(RF(SMax-k)))*
(1 - Eȓ(SMax-k+1,α)+[VR(SMax-k+1, RF(SMax-k+1))])}

with optimality achieved at ᾶ=αR(SMax-k, RF(SMax-k)).

19

Appendix E. Derivation of Ruin Factor Bucket Probabilities
Assume we are at time t and RF(t) > 0 is known.
P(RF(t+1) in Bucket #1)
=

P(0 < RF(t+1) ≤ (1.5)*(1/PR))

=

P(0 <

=

P( >

=

P ∞>ȓ

=

P(RF t ∗ 1

=

1 - Fȓ(t+1,α)(RF t ∗ 1

ȓ
ȓ

(E.1a)

≤ (1.5)*(1/PR) )

,
,

)

. ∗

,

RF t ∗ 1

(E.1c)
)

(E.1d)

< ∞)

(E.1e)

. ∗

≤ ȓ

. ∗

(E.1b)

,

)

. ∗

(E.1f)

P(RF(t+1) in Bucket #i), for (2 ≤ i ≤ (PR)*RFMax)
=

P(i/PR - 1/(2PR) < RF(t+1)) ≤ i/PR + 1/(2PR))

(E.2a)

=

P((1/PR)*(i - 1/2) <

(E.2b)

=

P(RF(t)*(1

ȓ

≤ (1/PR)*(i + 1/2))

,

>ȓ

∗

(E.2c)

,

RF(t)*(1
=

P(RF(t)*(1

ȓ

∗

(E.2d)

,

RF(t)*(1
=

Fȓ(t+1,α)(RF(t)*(1

∗

)

∗

–

Fȓ(t+1,α)(RF(t)*(1

20

)

∗

∗

)

(E.2e)

P(RF(t+1) in Bucket #(PR)*RFMax+1)
=

P(RFMax + 1/(2PR) < RF(t+1) < ∞)

=

P(RFMax + 1/(2PR) <

=

P(

=

P(RF t ∗ 1

=

P(RF t

=

Fȓ(t+1,α)(RF t ∗ 1

ȓ

,

ȓ

,

< ∞

(E.3b)

)

(E.3c)

/

ȓ

ȓ

)

/

(E.3a)

,

> RF t )

(E.3d)

))

(E.3e)

RF t ∗ 1

,

/

)) – Fȓ(t+1,α)(RF(t))

/

(E.3f)

Lastly, in (14) these probabilities are conditional on RF(t+1) > 0 ↔ ȓ(t+1,α) > RF(t), therefore
each must be divided by P(ȓ(t+1,α) > RF(t)) = 1 - Fȓ(t+1,α)(RF t

for their sum to equal 1 and

constitute a valid PMF. What we have derived is the conditional probability that the next ruin
factor is in any of the pre-defined buckets, given that no ruin occurs at the next time point. To
summarize the above, at time t the probability that RF(t+1) falls in bucket i, i = 1, 2, …,
(PR)*RFMax, (PR)*RFMax+1 given that RF(t+1) > 0 is defined by the PMF:
P(RF(t+1) in Bucket i | ȓ(t+1,α) > RF(t)) =
ȓ

,

∗
ȓ

ȓ

,

.

∗

,

,
∗

–

∗

ȓ
ȓ

,

∗

i=1

ȓ

,

∗

∗

,

/

ȓ

(E.4)

–

ȓ

,

,

,

21

,

i = 2, 3, …, (PR)*RFMax

i = (PR)*RFMax + 1

Appendix F. Derivation of Historical Stock and Bond Distributions
The hypothesis that real total stock (S&P 500) and bond (10 year T-Bond) returns
originate from normal distributions cannot be rejected (Anderson-Darling p-values of 0.707 and
0.243, respectively). Other distributions also provide an acceptable fit, such as the lognormal,
however only using a non-zero location parameter in the case of real stock returns. We thus
assume that the underlying distributions for real total returns at year t are given by, N(µ, σ):
S&P 500 Returns = r(s,t) ~ N(0.0825, 0.2007)

(F.1a)

10 year T-Bond Returns = r(b,t) ~ N(0.0214, 0.0834)

(F.1b)

In addition, a small positive correlation of ρ = 0.04387 was measured between real stock
and bond returns at the same time point, which will be carried through the analysis. We note that
these are unconditional distributions and would roughly trace out the shape of corresponding
histograms plotted using these returns. If the returns from either investment exhibit serial
correlation they cannot be treated as random samples from their unconditional distributions. It is
well known, for example, that the inflation rate exhibits strong serial correlation.
To determine the nature of serial correlation within real stock and bond returns, we
assume that the average return for each is stable, or that the processes are stationary. We
estimate the overall mean for each process by the sample average, which is subtracted from the
returns leaving them as deviations from their respective means. A stationary process that is
autoregressive of order p, denoted AR(p), is one that explicitly models the current observation as
a linear function of the past p observations, under the assumption that they are correlated. The
general AR(p) model is defined as (Box et al. (1994)),
Yt = ϕ1Yt-1 + ϕ2Yt-2 + … + ϕpYt-p + at
where Yt reflects the centered observation at time t and at ~ N(0, σa2) are iid error terms.

22

(F.2)

Note that we assume the Y’s are correlated despite the fact that the error terms are iid.
This is because, for instance, Yt and Yt+1 are both functions of the same random error term at. To
determine the autoregressive order of a stationary process we examine the autocorrelations and
partial autocorrelations. The autocorrelations at lag k are estimated by the sample correlations of
all observations k time points apart and the partial autocorrelations at lag k reflect the correlation
of observations k time points apart, after accounting for the correlation at lags < k. If an AR(k)
model is fit to the data, the estimate of ϕk is also an estimate of the k-th partial autocorrelation. It
is appropriate to examine both of these statistics up to lag N/4, where N reflects the number of
data points collected (Box et al. (1994), here N=86). Two noteworthy findings are that both the
autocorrelations and partial autocorrelations at lags > p are approximately N(0,

1/N) in an

AR(p) process for large N (Box et al. (1994)). Since approximately 95% of the data in a normal
distribution lies within 2 standard deviations of the mean, this result can be used to assess the lag
at which serial correlation ends in a stationary process (Box et al. (1994)).
The autocorrelation estimates for inflation-adjusted stock returns are shown below on the
left side of Figure F1.

None of the 21 autocorrelations exceed the 2 1/N (= 0.21567)

threshold, represented by the dashed line. This indicates that an AR(0) model is appropriate.
The partial autocorrelations were estimated by fitting successive AR(k) models, for k=1, 2, …,
21, where the ϕ-estimates for each model were taken as solutions to the Yule-Walker equations
(Box et al. (1994)). The partial autocorrelation estimates are shown on the right side of Figure
F1 and confirm the finding that real stock returns do not exhibit serial correlation. Since the
autocorrelation estimates shown (not partial) and their large sample distribution apply to any
linear stationary process, including a moving average MA(q), or mixed autoregressive-moving
average ARMA(p, q) process, we conclude that these returns originate as completely random
23

samples. Note that autocorrelations tail off after p in an AR(p) process, but drop off abruptly
after q in an MA(q) process (Box et al. (1994)).
Figure F1. Analysis of Serial Correlation in Real Stock Returns

We perform a similar analysis on real bond returns and display the results in Figure F2
below.

None of the 21 autocorrelations or partial autocorrelations exceed the 2 1/N (=

0.21567) tolerance limits indicating that an AR(0) process is appropriate for these returns.
Again, since the autocorrelation estimates and their large sample distribution apply to any linear
stationary process, we conclude that these returns also originate as completely random samples.
Figure F2. Analysis of Serial Correlation in Real Bond Returns

24

We note that these findings coincide with the theory that markets are efficient (at least in
the weak sense), meaning that prices adjust for any predictive capacity inherent in historical
patterns, thereby instantly removing the predictive information. Lastly, we assume that real total
stock and bond returns jointly follow a bivariate normal distribution.
We can test the stationarity assumption after fitting an AR(p) model to the data. The
fitted model defines a characteristic polynomial in the backshift operator B (where Bk[Yt]=Yt-k)
using the estimated coefficients ϕ and satisfying ϕ(B)*Yt = at from (F.2). We immediately see
that ϕ(B) = 1 – ϕ1B – ϕ2B2 – … – ϕpBp is such a polynomial and stationarity is confirmed if all
roots have magnitude > 1 (Box et al. (1994)). (Complex roots would appear as conjugate pairs, a
± bi, with magnitude √a

b .) Since an AR(0) process was found appropriate for both

returns, the stationarity test was not performed. Note that MA(q) models are stationary by
construct, and ARMA(p, q) models are stationary when their autoregressive component meets
the condition just stated (Box et al. (1994)). Regression models were also fit and indicate that
the average real return for stocks and bonds does not change as a linear function of time.
Due to the nature of compounding, investors generally do not experience the arithmetic
mean return of a security over time. The geometric mean is a more accurate reflection of the
compounded average performance. Based on the 86 years of historical data used here, the
geometric means of inflation-adjusted stock and bond returns are 0.063 and 0.018, respectively.
Eighty-six years of data were simulated (N=1 million times each) using the distributions for
stocks and bonds given in (F.1a) and (F.1b) and the average geometric mean values were 0.063
and 0.018, respectively, matching their historical counterparts.6

6

Both negative and positive outliers were removed in the same relative proportions to avoid square roots of
negative numbers. (About 1 outlier per 70,000 simulated histories was removed.)

25

At each stage/state of the discrete DP (see Figure 5) the minimum probability of ruin is
found over all asset allocations in the set α{·}. Let R(s,t) and R(b,t) denote the actual total returns
(not inflation-adjusted) at year t, for stocks and bonds respectively. A portfolio with stock and
bond proportions α, and (1-α) thus has a total return represented by the RV,
R(t,α) = α*R(s,t) + (1-α)*R(b,t),

(F.3a)

1 + R(t,α) = 1 + [α*R(s,t) + (1-α)*R(b,t)].

(F.3b)

with compounded return:

Adjusting for inflation, the compounded real return for this portfolio at year t is given by the RV:
1 + r(t,α) | It, R(t,α) = [1 + R(t,α)]/(1+It) = {1 + [α*R(s,t) + (1-α)*R(b,t)]}/(1+It)

(F.4a)

→ 1 + r(t,α) | It, R(t,α) = {α*[1 + R(s,t)] + (1-α)*[1 + R(b,t)]}/(1+It)

(F.4b)

→ 1 + r(t,α) | It, R(t,α) = α*[1 + R(s,t)]/(1+It) + (1-α)*[1 + R(b,t)]/(1+It)

(F.4c)

→ 1 + r(t,α) | It, R(t,α) = α*(1 + r(s,t)) + (1-α)*(1 + r(b,t))

(F.4d)

→ 1 + r(t,α) | It, R(t,α) = α*N(1.0825, 0.2007) + (1-α)*N(1.0214, 0.0834)

(F.4e)

→ 1 + r(t,α) | It, R(t,α) ~ N(α*(1.0825) + (1-α)*(1.0214),

(F.4f)

α ∗ 0.2007

1

α

∗ 0.0834

2∗α∗ 1

α ∗ .00073 ),

where in general,
Var(αX1 + (1- α)X2) = α2*Var(X1) + (1-α)2*Var(X2) + 2*α*(1- α)*Cov(X1,X2)

(F.4g)

and,
Cov(X1,X2) = ρ* Var X

∗ Var X ,

(F.4h)

holds for any 2 random variables X1 and X2. Finally, assuming an expense ratio of ER is paid to
the financial institution yearly, the inflation/expense-adjusted return at year t is given by the RV:
ȓ(t,α) = (1+r(t,α))(1-ER) ~ N((1-ER)*[α*(1.0825) + (1-α)*(1.0214)],
(1-ER ∗ α ∗ 0.2007

1

α
26

∗ 0.0834

2∗α∗ 1

(F.5)
α ∗ .00073 ) .

Appendix G. Miscellanea
G.1 The Equity Glide-Path Debate
Much of the pracitioner-based literature on retirement research has been devoted to
finding the best equity glide-path to use in decumulation. We have summarized the conclusions
from several such studies within the literature review section of this paper. In fact a spirited
debate on the topic has ensued and a natural question to ask is whether or not our model can
ultimately settle the issue. Let αt ∈ [0, 1] denote the equity ratio used at time t for t=0, 1, …, TD-1
(fixed). The inflation/expense-adjusted random return between times t=0 and t=1 is ȓ

which

,

has the asset allocation set at time t=0 and is observed at time t=1. We define an equity glidepath as any set α{gp} = {α0, α1, …, αTD-1} which is fixed at time t=0. Common equity glide-paths
would then satisfy the conditions: declining (αi > αj), rising (αi < αj), and constant (αi = αj) ∀ i < j =
0, 1, …, TD-1.

Applying our notation, the probability of ruin at any time point during

decumulation using equity glide-path α{gp} is:
P(Ruin(> 0))
= P(Ruin(1) ∪ Ruin(2) ∪ … ∪ Ruin(TD))

(G.1a)

= 1 – P(RuinC(1) ∩ RuinC(2) ∩ … ∩ RuinC(TD))

(G.1b)

= 1 – P( ȓ

,

> RF(0)] ∩ ȓ

> RF(1)] ∩ … ∩ ȓ

,

> RF(TD-1)])

,

(G.1c)

Independent RVs for any equity glide-path α{gp}.
= 1–

…

ȓ

,

,ȓ

,

,…,ȓ

,

ȓ

,

(G.1d)
ȓ

,

…

ȓ

,

The joint PDF of all future returns
using equity glide-path α{gp}.
The expression in (G.1d) is the probability of ruin at any time point during decumulation given
equity glide-path α{gp}. For any given α{gp} the code that implements our models (see Appendix
27

H) can be used to derive this probability with minimal modification (by fixing αt to its known
value at each time point). Since we are evaluating a given strategy, simulation can also be used.
To enter the equity glide-path debate we turn this problem around and seek the set of values {α0,
α1, …, αTD-1} that minimizes the expression in (G.1d). We recognize this as the optimization
problem addressed in Section II-G but with the added constraint α(t, RF1(t)) = α(t, RF2(t)) ∀ t = 0,
1, …, TD-1 and 0 < RF1(t) < RF2(t). This constraint ensures that the asset allocation αt at time t is
constant. Proceeding by induction as before the step at time t=TD does not change, but we
encounter a problem at time t=TD-1. Assuming RF(TD-1) is known there is no single αTD-1 that
minimizes P(Ruin(TD)) ∀ RF(TD-1) > 0. Without imposing new conditions or rules we cannot
make the decision at time t=TD-1 and therefore cannot proceed to time t=TD-2.
While our model cannot help settle the equity glide-path debate, we see that framing the
retirement decumulation problem this way is flawed. The retiree is under no obligation to fix
α{gp} at time t=0 and the models we have proposed in (9) and (12) do not. We delay the asset
allocation decision until each time point t when RF(t) is known and express the optimal α to use
between times t and t+1 as a function of (t, RF(t)), namely ᾶ=α(t, RF(t)). We have also shown
that the ruin factor, RF(t), encapsulates into a single value all prior information needed to make
this decision. Further, since adding constraints to an optimization problem cannot improve the
objective, our model will perform at least as well as any based on a predetermined equity glidepath α{gp}. To achieve optimality with a predetermined α{gp} one must somehow guess in
advance which glide-path will be right. Simply assessing all of them is not sufficient. (Note that
with TD=30 and Pα=1,000 there are 1,00030 possible equity glide-paths.) We conclude that all
solutions based on preset equity glide-paths, including those used in target-date funds, are
suboptimal when implemented in decumulation with a safe withdrawal rate WR.

28

G.2 Splitting the Joint PDF during Backward Induction
In Appendix G.1 we highlighted a difference between this research and much of the
literature. Namely, we delay the asset allocation decision until the last moment, which allows us
to compute and use the ruin factor RF(t) when making the decision. The optimal α at each time
point is denoted by ᾶ=α(t, RF(t)) and our goal is to determine this function for each t and RF(t) > 0.
When the retiree sets their asset allocation at time t the observed return at time t+1 is assumed
independent of all returns observed at times ≤ t. During backward induction, however, we split
the PDF of the time t+1 unobserved return from the joint PDF of all time > t+1 unobserved
returns. At time t we compute RF(t) then optimize to find ᾶ=α(t, RF(t)). This asset allocation
defines the PDF of the RV ȓ(t+1,ᾶ) which is the inflation/expense-adjusted return for the next time
point. The subsequent time point’s return is ȓ(t+2,ᾶ) where ᾶ=α(t+1, RF(t+1)), but since RF(t+1) =
RF(t)/[ȓ(t+1,ᾶ) – RF(t)], the PDF of the RV ȓ(t+2,ᾶ) is a function of the RV ȓ(t+1,ᾶ), raising the question:
Can the joint PDF be separated? We can separate the joint PDF by conditioning on ȓ(t+1,ᾶ) and we
know the marginal PDF of ȓ(t+1,ᾶ) from independence. This will be made clear by first deriving
the value function using a different approach and then returning to our original method and
demonstrating precisely the conditioning that occurs, allowing separation of the joint PDF.
Without loss of generality, we will focus on the induction from Appendix C.4, where TD is fixed
and t=TD-3. We observe that for any event and random variable X ~ fX(x) (see, Ross (2007)):
P(Event) =

P Event | X

∗

d

(G.2a)

During induction at time t=TD-3 from Appendix C.4 we are interested in:
P(Ruin(> TD-3)) = P(Ruin(TD-2) ∪ Ruin(TD-1) ∪ Ruin(TD))

(G.2b)

= 1 - P(RuinC(TD-2) ∩ RuinC(TD-1) ∩ RuinC(TD))

(G.2c)

= 1 - P(RuinC(TD-2))*P(RuinC(TD-1) ∩ RuinC(TD) | RuinC(TD-2))

(G.2d)

29

The value function is:
1 - P(RuinC(TD-2))*

V(TD-3, RF(TD-3)) = Min

C

C

(G.2e)

C

P(Ruin (TD-1) ∩ Ruin (TD) | Ruin (TD-2))
Applying the definition of
a conditional probability.

→ V(TD-3, RF(TD-3))

1 - P(RuinC(TD-2))*

= Min

∩

(G.2f)

∩

Since RuinC(TD-2) ≡ ȓ(TD-2,α) > RF(TD-3).

→ V(TD-3, RF(TD-3))

(G.2G)
= Min

1 - P(ȓ(TD-2,α) > RF(TD-3))*
∩

∩ ȓ TD 2,α

ȓ TD 2,α

RF TD 3

RF TD 3

We now apply the observation in (G.2a) to the numerator in (G.2G), which yields:
P(RuinC(TD-1) ∩ RuinC(TD) ∩ ȓ(TD-2,α) > RF(TD-3))
=

P Ruin T

1 ∩ Ruin T

∩ ȓ TD

2, α

(G.2H)

RF TD

3 | ȓ TD

2, α

∗

ȓ TD 2,α

d

We separate this integral into 2 components, where the first vanishes and the second simplifies:
This probability is zero when ȓ(TD-2,α) ≤ RF(TD-3).
(G.2I)
=

P Ruin T

1 ∩ Ruin T

∩ ȓ TD 2, α

RF TD

3 | ȓ TD

2, α

∗

+

P Ruin T

1 ∩ Ruin T

∩ ȓ TD 2, α

RF TD

3 | ȓ TD

2, α

∗

ȓ TD 2,α
ȓ TD 2,α

d
d

This event will always occur when ȓ(TD-2,α) > RF(TD-3) and can be dropped.
In (G.2G), (G.2I) was the numerator of an expression [·], and including the denominator yields:
[·] =

∩

| ȓ TD 2,α
ȓ TD 2,α

30

∗

ȓ TD 2,α

(G.2J)

This is the PDF of the RV ȓ(TD-2,α)+ = (ȓ(TD-2,α) | ȓ(TD-2,α) > RF(TD-3)).
(G.2K)
=

ȓ TD 2,α
ȓ TD 2,α

∗ P Ruin T

1 ∩ Ruin T

| ȓ TD

2, α

= Eȓ(TD-2,α)+[ P(RuinC(TD-1) ∩ P(RuinC(TD) | ȓ(TD-2,α)]

d
(G.2L)

Substituting (G.2L) back into (G.2G) clarifies our task, which is to choose an α that minimizes
the bracketed expression {·} and assign it to the function V(TD-3, RF(TD-3)).
→ V(TD-3, RF(TD-3))
= Min

(G.2M)
1 - P(ȓ(TD-2,α) > RF(TD-3))*
Eȓ(TD-2,α)+[P(RuinC(TD-1) ∩ P(RuinC(TD) | ȓ(TD-2,α)]

The events RuinC(TD-1) and RuinC(TD) are determined by future returns, namely those observed
at times t=TD-1 and t=TD. Since they are based on future asset allocation decisions we are free to
consider any values. However, to minimize the expression we must use the largest possible
value of [P(RuinC(TD-1) ∩ P(RuinC(TD) | ȓ(TD-2,α)] for each α. If we do not, the bracketed term
{·} cannot be a minimum. We recognize that [P(RuinC(TD-1) ∩ P(RuinC(TD) | ȓ(TD-2,α)] is
maximized at [1 – V(TD-2, RF(TD-2))] by definition of the value function and (G.2M) becomes:
→ V(TD-3, RF(TD-3))
= Min

(G.2N)
1 - (1 – Fȓ(TD-2,α)(RF(TD-3)))*
(1 - Eȓ(TD-2,α) +[V(TD-2, RF(TD-2))])

This is identical to what was derived in (C.9). In (C.8c), during induction for fixed TD at time
t=TD-3, we express a probability statement as the integration over a corresponding joint PDF
f(ȓ(TD-2,α), ȓ(TD-1,ᾶ), ȓ(TD,ᾶ)), then split f(ȓ(TD-2,α)) from f(ȓ(TD-1,ᾶ), ȓ(TD,ᾶ)) in (C.8e). Since ȓ(TD-1,ᾶ) uses
ᾶ=α(TD-2, RF(TD-2)) and RF(TD-2) = RF(TD-3)/[ȓ(TD-2,α) - RF(TD-3)] some may question the
validity of this operation. The same is true for RVs ȓ(TD-2,α) and ȓ(TD,α). The numerator of the ratio
31

in (C.8c) is P[ȓ(TD-1,ᾶ) > RF(TD-2) ∩ ȓ(TD,ᾶ) > RF(TD-1) ∩ ȓ(TD-2,α) > RF(TD-3)]. Using (G.2a) we can
express this probability as:
P[ȓ(TD-1,ᾶ) > RF(TD-2) ∩ ȓ(TD,ᾶ) > RF(TD-1) ∩ ȓ(TD-2,α) > RF(TD-3)]
=

(G.2O)

P(ȓ(TD-1,ᾶ) > RF(TD-2) ∩ ȓ(TD,ᾶ) > RF(TD-1) ∩ ȓ(TD-2,α) > RF(TD-3) | ȓ(TD-2,α) = r)*

d

ȓ TD 2,α

This probability is zero when ȓ(TD-2,α) ≤ RF(TD-3).
(G.2P)
=

P(ȓ(TD-1,ᾶ) > RF(TD-2) ∩ ȓ(TD,ᾶ) > RF(TD-1) ∩ ȓ(TD-2,α) > RF(TD-3) | ȓ(TD-2,α) = r)*

ȓ TD 2,α

d

+

P(ȓ(TD-1,ᾶ) > RF(TD-2) ∩ ȓ(TD,ᾶ) > RF(TD-1) ∩ ȓ(TD-2,α) > RF(TD-3) | ȓ(TD-2,α) = r)*

ȓ TD 2,α

d

This event will always occur when ȓ(TD-2,α) > RF(TD-3) and can be dropped.
=

∗P(ȓ(TD-1,ᾶ) > RF(TD-2) ∩ ȓ(TD,ᾶ) > RF(TD-1) | ȓ(TD-2,α) = r) d

ȓ TD 2,α

(G.2Q)

This probability is replaced by its integral expression.
(G.2R)
=

ȓ

,

ȓ

,ᾶ

,ȓ

,ᾶ

ȓ

,ᾶ

ȓ

,ᾶ

ȓ

,

The joint PDF has been separated which
validates the operation in (C.8e).
We conclude by emphasizing some subtle points from above. At time t we set the α that
defines the RV ȓ(t+1) which is observed at time t+1. At time t+1 we set the α that defines the RV
ȓ(t+2) which is observed at time t+2. Regardless of the ȓ(t+2) we choose to represent our portfolio
between times t+1 and t+2 it will be independent of ȓ(t+1) via market efficiency. However, the
choice of ȓ(t+2) does depend on ȓ(t+1). Further, since we observe returns sequentially in time we
can always condition on knowledge of what has happened previously if it helps to simplify
expressions. This is true because backward induction assumes the prior RF(t) is known, and it is
also true when we use the model since we start at time t=0 and proceed sequentially in time.
32

Appendix H. Full C++ Implementation
The implementation consists of 5 functions and one header that are compiled into a 64-bit
console application (executable). The application is launched by either double-clicking the
executable or by invoking it from the command line or via a batch file. A single parameter is
optionally accepted which is the number of tasks the user wants to process concurrently as it
runs. If no value is supplied the code determines the maximum number of concurrent processing
units on the machine running it and uses this value. The maximum number of concurrent
processing units is generally the number of CPU cores. If an Intel® processor with HyperThreading® (HT) is used the capacity of each core doubles. For example, the runtime statistics
given next are from a Dell® Precision® workstation with dual 10-core 2.5GHz Intel® Xeon®
processors. Since all Xeon® processors use HT this PC can execute 2x10x2=40 concurrent jobs.
The application is multi-threaded and the task of processing each time point is split as follows.
First the code approximates the RF(t) bucket at which heavy algorithm pruning begins. It then
splits the remaining buckets into equal-sized collections and launches jobs to process all buckets
simultaneously in separate threads. Once all threads for a given time point finish the results are
aggregated and processing for the next time point begins.
When launched the application queries the user for the location of a directory that
contains a control file and optionally a file of age probabilities (required for random TD). The
control file defaults to filename “control.txt” and the age probabilities file to “ageprobs.txt”.
Both filenames can be changed in the header if desired. Examples of 6 control files (including
runtime statistics) and the age probability file (derived from SSA.gov) are shown next. This
application uses the C++ boost© library which will need to be installed (including binaries) on
the machine creating the executable and it was optimized for speed when compiled using
Microsoft® Visual Studio® (a project level setting).
33

Input File: control.txt
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Contents: (A 3-line text file of space-delimited values as shown below.)
StockMean StockVar BondMean BondVar StockBondCov RFMax ER PrunePwr
PR Pα
NumRand Details

-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Example 1: (Fixed TD=50, PR=1,000, Pα=100, RFMax=2.75, ER=0.0%, Runtime=00:00:43)
0.082509 0.0402696529 0.021409 0.0069605649 0.0007344180 2.75 0.000 4.00
1000 100
0 50

-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Example 2: (Fixed TD=50, PR=5,000, Pα=1,000, RFMax=2.75, ER=0.5%, Runtime=00:28:02)
0.082509 0.0402696529 0.021409 0.0069605649 0.0007344180 2.75 0.005 4.00
5000 1000
0 50

-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Example 3: (Fixed TD=50, PR=10,000, Pα=1,000, RFMax=2.75, ER=1.0%, Runtime=01:42:36)
0.082509 0.0402696529 0.021409 0.0069605649 0.0007344180 2.75 0.010 4.00
10000 1000
0 50

-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Example 4: (Random TD/N=2, SMax=46, PR=1,000, Pα=100, RFMax=2.75, ER=0.0%, Runtime=00:01:17)
0.082509 0.0402696529 0.021409 0.0069605649 0.0007344180 2.75 0.000 4.00
1000 100
2 M 65 F 67

-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Example 5: (Random TD/N=9, SMax=50, PR=5,000, Pα=1,000, RFMax=2.75, ER=0.5%, Runtime=02:22:08)
0.082509 0.0402696529 0.021409 0.0069605649 0.0007344180 2.75 0.005 4.00
5000 1000
9 M 62 F 63 M 66 F 67 F 70 F 72 M 74 M 75 F 76

-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Example 6: (Random TD/N=1, SMax=48, PR=10,000, Pα=1,000, RFMax=2.75, ER=1.0%, Runtime=11:01:24)
0.082509 0.0402696529 0.021409 0.0069605649 0.0007344180 2.75 0.010 4.00
10000 1000
1 F 65

-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Note: Runtime units are HH:MM:SS. The means/variances match the historical values that
were derived in Appendix F. A 50-year fixed TD solution contains all solutions for fixed TD
under 50 years. The parameter PrunePwr=4.0 indicates that heavy algorithm pruning begins
when P(Ruin) equals the largest possible P(Ruin) for that time point to within 4 decimal places.
Once ruin becomes a near certainty α=1 and it is not necessary to evaluate all α-values. A
higher PrunePwr value has a longer runtime and the user can confirm an equivalent solution by
adjusting this value. The # years processed for random TD is the larger of (111-min male age)
and (113-min female age) where 111 and 113 are the maximum possible male/female ages,
respectively, as derived from SSA.gov life-tables. (See the file “ageprobs.txt” below.) When TD
is random the last possible withdrawal is at t=SMax but the last decision point is at t=SMax-1.
34

Input File: ageprobs.txt
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

50 0.005346265174752670 0.003283385442263650
51 0.005810688210135220 0.003544801480660440
52 0.006264310709811210 0.003795760877521360
53 0.006717933209487190 0.004015350349774660
54 0.007171555709163170 0.004234939822027960
55 0.007657579815958870 0.004475442577353010
56 0.008176005529874280 0.004747315257285670
57 0.008683630708083120 0.005071471144897680
58 0.009180455350585390 0.005458366881724930
59 0.009677279993087660 0.005897545826231530
60 0.010228107314122800 0.006399464619953360
61 0.010822136777984200 0.006943209979818680
62 0.011513371063204700 0.007518325264291610
63 0.012269408562664700 0.008103897190300420
64 0.013111850347777200 0.008741752323988580
65 0.014029895882835800 0.009452803948427840
66 0.015034345703546900 0.010247508705154100
67 0.016038795524258000 0.011104953311095500
68 0.017064846416382300 0.012035594407788100
69 0.018112498379919600 0.013049888636767600
70 0.019278956236229300 0.014168749281105900
71 0.020531818378191600 0.015392176340802900
72 0.021827882662980100 0.016667886608179200
73 0.023134747483475200 0.017943596875555500
74 0.024463213375383400 0.019282046992147100
75 0.025899684624357400 0.020787803373312600
76 0.027422560158983900 0.022429496094444400
77 0.028891433015077500 0.024050275532504500
78 0.030273901585518600 0.025629228404421100
79 0.031591566941720300 0.027208181276337700
80 0.032887631226508800 0.028891700563613000
81 0.034140493368471100 0.030721612832390500
82 0.035166544260595300 0.032561981742703900
83 0.035922581760055300 0.034339610803802000
84 0.036365403724024700 0.036033586732613200
85 0.036516611223916700 0.037539343113778700
86 0.036300600509785300 0.038804596739619200
87 0.035663368903097600 0.039693411270168200
88 0.034572514796733900 0.040132590214674900
89 0.032963234976454800 0.040028023799316100
90 0.030857130513673500 0.039285602250269300
91 0.028265001944096400 0.037894868925998300
92 0.025305655160496000 0.035845367184967500
93 0.022076294984231200 0.033189380234856200
94 0.018717328379487600 0.030021017849487100
95 0.015347561239037500 0.026371649953467900
96 0.012129001598479300 0.022481779302123700
97 0.009277660171944530 0.018591908650779500
98 0.006825938566552900 0.014890257547081000
99 0.004881842139370110 0.011585958821745600
100 0.003391368211863310 0.008804492173203810
101 0.002300514105499630 0.006545857601455560
102 0.001522875534626520 0.004726401974213920
103 0.000972048213591394 0.003335668649943010
104 0.000604829999567979 0.002269091213284120
105 0.000356417678316845 0.001495299739629630
106 0.000216010714131421 0.000951554379764307
107 0.000108005357065710 0.000575115284472933
108 0.000064803214239426 0.000345069170683760
109 0.000032401607119713 0.000188219547645687
110 0.000010800535706571 0.000094109773822844
111 0.000010800535706571 0.000052283207679358
112 0.000000000000000000 0.000020913283071743
113 0.000000000000000000 0.000010456641535872

35

Note: Columns in this file represent:
age, male death probability, and female
death probability for a 50-year old.
They are computed from life-tables
published at SSA.gov and each column
sums to 1. For example, a male aged
50 has a 0.005346265174752670
probability of death before age 51. As
seen, the maximum possible male age is
111 and the maximum possible female
age is 113. There are no blank lines in
this file and each field is delimited by a
single space.

Code File: minpruin.cpp
/*
/ Copyright (c) 2014 Chris Rook, All Rights Reserved.
/
/
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
/
/
1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
/
2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the
/
documentation and/or other materials provided with the distribution.
/
3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this
/
software without specific prior written permission.
/
/
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
/
TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
/
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
/
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
/
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
/
EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. (Source: http://opensource.org/licenses/BSD-3-Clause.)
/
/ Filename: minpruin.cpp (Defines the entry point for the console application.)
/
/ Function: main()
/
/ Summary:
/
/
Compiling the code will generate a single executable file that can be invoked with a double click or via a batch file with one parameter. When invoked
/
the directory location where input files reside and output files are written to is queried, stored as a string vbl and passed to other functions.
/
The control file (control.txt) is then read. If TD is specified as random the function derhrates() is invoked which reads the file of male/female
/
death probabilities (ageprobs.txt) and derives the hazard rates that apply to the group listed in the last line of the control file. The hazard rates
/
are stored locally and written to the file specified by the constant hrfile for reference. When TD is fixed, hazard rates are not needed or derived.
/
When invoked, the executable may have either zero or one argument(s). (See below.) The arrangement specified by the control file determines the
/
total number of years to process. Iteration always begins at the terminal year (i.e., total # of years minus 1) and ends at year 0. This reflects the
/
fact that the first asset allocation (alpha) decision point is at time t=0 and the last is at time t=(total # years - 1). The first withdrawal attempt
/
(assuming survival) is at time t=1 and the last withdrawal attempt (assuming survival) is at time t=(total # years). For example, if TD is fixed
/
and 50 years are specified in the control file then processing starts at year=49 and ends at year=0. (This results in 50 alpha decisions made by the
/
retiree, the first at year=0 (i.e., t=0) and the last at year=49 (i.e., t=49). The total # of years to process for random TD depends on the specific
/
ages listed for the group (MPU). For example, if random TD for a 65 year old male and female couple is specified in the control file then a total of 48
/
years are processed since the maximum age allowed for a male is 111 and for a female is 113 (see the file ageprobs.txt). Processing 48 years of data
/
implies that 48 alpha decisions are made, here the first is at year=0 (both M/F at 65) and the last at year=47 (F is age 112). (See the paper for a
/
warning about dated hazard risk. Rerun the optimization if the hazard probabilities differ notably over time from those written to hrates.txt.) Using a
/
different ageprobs.txt file can change total # of years that need to be processed which, in general is: Max((Max male age - Min male age),(Max female
/
age - Min female age)). Here maximums are allowable and minimums are actual as specified in the control file. When processing each year the function
/
getprprobs() is invoked and reads the probabilites derived for the previous year. These are then passed to the function optimize() which derives the
/
probabilities and alphas for the current year and writes the results to a file for reference. Once complete, the function combfiles() is invoked which
/
combines all separate files into a single file to be read in as the prior probabilities during processing of the next year/time point. When at year=0,
/
combfiles() also concatenates all data from all prior years into both a vertical text file and 2 horizontal csv files, where the csv files can be read
/
into a spreadsheet. This is helpful when simulating the optimal solution to verify the result and only the csv file of optimal alphas is needed.
/
/
Notes: 1.) The terms "year" and "time point" are used interchangably throughout this application and mean the same thing.
/
2.) The file "ageprobs.txt" was derived from life tables at SSA.gov.
/
/ Parameters:
/
/
The executable takes zero or one argument(s) when invoked. If an argument is specified it reflects the maximum number of jobs to process concurrently
/
by the application as it executes. If no argument is specified this value defaults to the number of independent processing units that exist on the
/
machine running the application. (This application was developed and tested on the Windows 7® operating system.)
/
/ Input Files:
/
/
1.) Parameter control file of specific form whose name is defined by the string constant paramfile (defined in the header and located in the directory

36

/
specified by the user).
/
2.) Age probability file of specific form whose name is set as the string constant ageprfile (defined in the header file and located in the directory
/
specified by the user. This file is read and processed by the function derhrates()).
/
/ Output Files:
/
/
1.) The derived hazard rates for the random arrangement specified in the control file. The hazard rate file is written by function derhrates() using
/
filename specified by string constant hrfile defined in the header to the directory specified by the user.
/
2.) FinalResults_V.txt (Vertical text file of all results from all years processed, written by combfiles() after year=0 finishes.)
/
3.) FinalAlphaResults_H.csv (Horizontal csv file of optimal alphas, written by combfiles() after year=0 finishes. Use this file to simulate the result.)
/
4.) FinalProbResults_H.csv (Horizontal csv file of optimal probabilities, written by combfiles() after year=0 finishes.)
/
/
Note: A number of temporary data files are written to the directory specified by the user during the processing of each year and across years. These
/
are removed by the application as it runs. They are created to avoid carrying large amounts of data inside arrays that are not needed. This was
/
a design decision made to speed up runtime. When processing a given year/time point, only the prior year's data are needed. Data from all other
/
years/time points is not used (by nature of DPs).
/
/ Return Value:
/
/
This function returns the integer value of 0 (success) or an error code (failure).
/------------------------------------------------------------------------------------------------------------------------------------------------------------------------*/
#include "stdafx.h"
int main(int argc, char *argv[])
{
// Declare variables.
int precisions[2], numrand, numyears, p=0, rc, pllproc=0;
long buckets[2], * * bktarys, bktsprun;
double params[8];
long double * hr, * prV;
string rootdir, trbin;
vector<long> unqbkts;
boost::thread * t;
// Get the directory location where setup files are stored, strip any leading/trailing blanks.
cout << "Enter the directory where the setup files reside (eg, c:\\mprsetup\\): " << endl;
cin >> rootdir; cin.get();
boost::algorithm::trim(rootdir);
cout << endl;
// Read in parameter control file. Load first line into array params[] and 2nd line into array precisions[].
// (3rd line is processed conditionally based on the arrangement specified by the user i.e., fixed vs random TD.)
ifstream getparams(rootdir+paramfile);
if (getparams.is_open())
{
getparams >> params[0] >> params[1] >> params[2] >> params[3] >> params[4] >> params[5] >> params[6] >> params[7];
getparams >> precisions[0] >> precisions[1];
getparams >> numrand;
if (numrand > 0)
{
char * genders = new char[numrand];
int * ages = new int[numrand];
if (numrand == 1)
getparams >> genders[0] >> ages[0];
// Gender and age of single person with random TD.
else
// Random TD for group. Load arrays with genders and ages.
{
string people;
std::getline(getparams, people);
stringstream ss(people);
for (int i=1; i <= numrand && ss.good(); ++i)
{

37

ss >> genders[i-1] >> ages[i-1];
genders[i-1] = toupper(genders[i-1]);
// Data check for invalid gender. (Invalid ages are checked in derhrates().)
if (genders[i-1] != 'M' && genders[i-1] != 'F')
{
cout << "ERROR: Invalid gender in file " + rootdir + paramfile + " for person #" << i << ": " << genders[i-1] << endl;
cout<< "EXITING...main()..." << endl; cin.get();
exit (EXIT_FAILURE);
}
}

}

// Build the hazard rates if random TD is specified and no file exists in the directory. (Don't leave an old one sitting there.)
ifstream gethrates(rootdir+hrfile);
numyears = derhrates(rootdir, numrand, genders, ages, gethrates.is_open());
hr = new long double[numyears+1];
delete [] genders; genders = nullptr;
delete [] ages; ages = nullptr;
// File of hazard rates has now been built and written for reference, load it to the hr[] array.
if (!gethrates.is_open())
gethrates.open(rootdir+hrfile);
// Read in survival probabilities that have been written to the file hrates.txt for this arrangement.
// (Doing it this way is useful when needing to run a single year for debugging purposes.)
if (gethrates.is_open())
{
while (!gethrates.eof() && p < numyears+1)
{
gethrates >> hr[p] >> trbin;
p=p+1;
}
gethrates.close();
}
else
{
cout << "ERROR: Could not read file: " << rootdir + hrfile << endl;
cout << "EXITING...main()..." << endl; cin.get();
exit (EXIT_FAILURE);
}
}
else
{

}
else
{

}

// For fixed TD load hr[] with zeros (hrates.txt is not created).
// Using a fixed TD is similar to using a random TD with all discrete hazard probabilities = 0 except the one at TD, which equals 1.
// That is, P(TD=t|TD>=t) = 0 for t=0,1,2,...,TD-1 and P(TD=t|TD>=t)=1 for t=TD. (Pr[death before TD] = 0 and Pr[death at TD]=1.)
getparams >> numyears;
// Account for TD using HR=1.00.
hr = new long double[numyears+1];
for (int j=0; j<numyears; ++j)
hr[j]=0.0;
hr[numyears]=1.00;

}
getparams.close();

cout << "ERROR: Could not read file: " << rootdir + paramfile << endl;
cout << "EXITING...main()..." << endl; cin.get();
exit (EXIT_FAILURE);

38

// Parse arguments to main(). Pllproc is optional, default is optimal # of threads on the machine running the application.
if (argc == 2)
pllproc = stoi(argv[1]);
else if (argc != 1)
{
cout << "ERROR: Parameter misspecification. Incorrect # of parameters to the executable (expecting zero or one...)." << endl;
cout << "EXITING...main()..." << endl; cin.get();
exit (EXIT_FAILURE);
}
// When pllproc==0, replace with the # of independent processing units on the computer running the application.
// Note that main() runs in its own thread but that is not accounted for here as it sits idle.
if (pllproc == 0)
pllproc = boost::thread::hardware_concurrency();
if (pllproc == 1)
pllproc = pllproc++;
// If just one processing unit use 2 threads.
// Iterate over all years, launching separate threads to proccess optimal size collections of buckets concurrently.
for (int yr=numyears-1; yr>=0; yr--)
{
cout << endl << "Processing for year " << yr << " has begun ..." << endl;
// Retrieve the prior year's probabilities which have been written to a file and load into an array.
prV = new long double[(long) (params[5]*precisions[0] + 1.5)];
unqbkts=getprprobs(rootdir, yr, numyears, params[5], precisions[0], prV, hr[yr+1]);
// Function call to determine pruning/parallel processing parameters for the current year.
buckets[0]=0;
buckets[1]=0;
rc = optimize(rootdir, yr, params, precisions, buckets, hr, prV, unqbkts);
// Derive the optimal # of buckets to process per run.
bktsprun = (long) (rc / (pllproc-1));
cout << "--> # Buckets processed per thread (excluding last one): " << bktsprun << endl;
// Build thread array of optimal size then launch calls to optimize() concurrently.
bktarys = new long * [pllproc];
t = new boost::thread[pllproc];
for (int i=1; i<=pllproc; ++i)
{
// Define buckets then process them in concurrent threads.
bktarys[i-1] = new long[2];
bktarys[i-1][0]=bktsprun*(i-1) + 1;
bktarys[i-1][1]=(i < pllproc)*bktsprun*(i) + (i==pllproc)*((long) (params[5]*precisions[0] + 0.5));
cout << "--> Begin concurrent processing of buckets " << bktarys[i-1][0] << " through " << bktarys[i-1][1] << " ..." << endl;
t[i-1] = boost::thread(optimize, rootdir, yr, boost::cref(params), boost::cref(precisions), boost::cref(bktarys[i-1]), boost::cref(hr),
boost::cref(prV), boost::cref(unqbkts));
}
// Wait for all threads to finish, then proceed.
for (int i=0; i<pllproc; ++i)
t[i].join();
// Free dynamic memory allocations.
for (int i=0; i<pllproc; ++i)
{
delete [] bktarys[i]; bktarys[i] = nullptr;
}
delete [] bktarys; bktarys = nullptr;
delete [] t; t = nullptr;
delete [] prV; prV = nullptr;

39

// Clear contents of unqbkts vector and reset capacity.
unqbkts.clear();
unqbkts.shrink_to_fit();

(This container is reused within the loop.)

// Combine all files for each year processed, and if at year=0 concatenate data files from all years.
combfiles(rootdir, yr, pllproc, bktsprun, params[5], precisions[0], numyears);
}

}

cout << "Processing for year " << yr << " has finished." << endl;

// Free memory of hazard rate array and exit.
delete [] hr; hr = nullptr;
return 0;

Code File: derhrates.cpp
/*
/ Copyright (c) 2014 Chris Rook, All Rights Reserved.
/
/
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
/
/
1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
/
2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the
/
documentation and/or other materials provided with the distribution.
/
3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this
/
software without specific prior written permission.
/
/
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
/
TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
/
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
/
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
/
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
/
EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. (Source: http://opensource.org/licenses/BSD-3-Clause.)
/
/ Filename: derhrates.cpp
/
/ Function: derhrates()
/
/ Summary:
/
/
Derive the discrete hazard rates for the random TD arrangement specified by the last line of the control file. This function is not called for
/
fixed TD. The derived hazard rates are written to a file whose name is specified by the string constant hrfile (defined in the header file), and
/
the total # of years that need to be processed are determined by this function and returned by the call.
/
/ Parameters:
/
/
1.) Root directory where input files reside and output files are written.
/
2.) # of persons specified by this arrangement in the control file.
/
3.) Array of genders for the random persons involved in this arrangement.
/
4.) Array of ages for the random persons involved in this arrangement.
/
5.) True/false indicator of whether or not the hazard rates file already exists. If it exists, this function only returns the total number of
/
years that need to be processed by this arrangement and the hazard rates are not re-derived. (Do not leave an old hazard rates file in the
/
root directory. If one is there it will be used but an alert message is given to the user. Kill the job with task manager if necessary.)
/
/ Input Files:
/
/
1.) Age probability file of specific form whose name/location is set as the string constant ageprfile which is defined in the header file and
/
resides in the directory specified by the user.
/
/ Output Files:

40

/
/
1.) The derived hazard rates for the random arragement specified in the control file using filename specified by hrfile which is defined in the
/
header file and resides in the directory specified by the user.
/
/ Return Value:
/
/
This function returns the # of years that need to be processed for this arrangement, representing the maximum # of alpha decisions to make.
/-------------------------------------------------------------------------------------------------------------------------------------------------------------*/
#include "stdafx.h"
int derhrates(const string root, const int npersons, const char * gndrs, const int * ags, const bool hrexst)
{
int age, prvage, remyrs, strtage=0, minmage=999, maxmage=0, maxamage, minfage=999, maxfage=0, maxafage;
long double prob1, prob2, * * perspmf, * * perscdf, sumprobs, mchksum=0, fchksum=0, * tdcdf, * hrates;
vector<long double> mprobs, fprobs;
// Read age probabilities for M/F into corresponding arrays. These are TD probabilities starting at a given age and will need to be adjusted for each person
// based on their age at retirement start. Also, use the file to derive the maximum possible M/F ages, which are ages that have non-zero probabilities.
ifstream getprobs(root+ageprfile);
if (getprobs.is_open())
{
while (!getprobs.eof())
{
getprobs >> age >> prob1 >> prob2;
if (strtage == 0)
{
strtage = age;
prvage = age-1;
}
if (age == prvage+1)
{
if (prob1 >= 0.0)
mprobs.push_back(prob1);
if (prob2 >= 0.0)
fprobs.push_back(prob2);
if (prob1 > 0.0)
maxamage=age;
if (prob2 > 0.0)
maxafage=age;
prvage = age;
}
}
getprobs.close();
}
else
{
cout << "ERROR: Could not read file: " << root + ageprfile << endl;
cout << "EXITING...derhrates()..." << endl; cin.get();
exit (EXIT_FAILURE);
}
// Confirm that the probabilities read sum to 1.00 for each gender (allow for floating pt precision diffs).
for (int i=(int) mprobs.size()-1; i>= 0; --i)
mchksum=mchksum+mprobs[i];
for (int i=(int) fprobs.size()-1; i>= 0; --i)
fchksum=fchksum+fprobs[i];
if (min(mchksum,fchksum) < (1.00 - pow(0.1,15)) || max(mchksum,fchksum) > (1.00 + pow(0.1,15)))
{
cout.setf(ios_base::fixed, ios_base::floatfield); cout.precision(25);
cout << "ERROR: Probabilities in " << root + ageprfile << " do not sum to 1 for one or both genders.
cout << "--> Male probability sum
= " << mchksum << endl;

41

See below..." << endl;

}

cout << "--> Female probability sum = " << fchksum << endl;
cout << "EXITING...derhrates()..." << endl; cin.get();
exit (EXIT_FAILURE);

// Find the minimum male/female ages specified in the control file at retirement start.
for (int i=0; i<npersons; ++i)
{
if (gndrs[i] == 'M' && ags[i] < minmage)
minmage = ags[i];
else if (gndrs[i] == 'F' && ags[i] < minfage)
minfage = ags[i];
}
// Find the maximum male/female ages specified in the control file at retirement start.
for (int i=0; i<npersons; ++i)
{
if (gndrs[i] == 'M' && ags[i] > maxmage)
maxmage = ags[i];
else if (gndrs[i] == 'F' && ags[i] > maxfage)
maxfage = ags[i];
}
// Check that
if (minmage <
{
cout
cout
cout
cout
cout
exit
}
// Check that
if (maxmage >
{
cout
cout
cout
cout
exit
}
// Check that
if (maxfage >
{
cout
cout
cout
cout
exit
}

no age specified on the last line of the control file is less than the start age in the probabilities file.
strtage || minfage < strtage)
<< "ERROR: An age in " << root + paramfile << " is less than the minimum age in " << root + ageprfile << "." << endl;
<< "--> Smallest male age (" << paramfile << ") = " << minmage << endl;
<< "--> Smallest female age (" << paramfile << ") = " << minfage << endl;
<< "--> Smallest allowable age (" << ageprfile << ") = " << strtage << endl;
<< "EXITING...derhrates()..." << endl; cin.get();
(EXIT_FAILURE);
no male age specified on the last line of the control file is greater than the last (>0 prob) male age in the age probs file.
maxamage)
<< "ERROR: A male age in " << root + paramfile << " is greater than the maximum allowed in " << root + ageprfile << "." << endl;
<< "--> Largest male age (" << paramfile << ") = " << maxmage << endl;
<< "--> Largest allowable male age (" << ageprfile << ") = " << maxamage << endl;
<< "EXITING...derhrates()..." << endl; cin.get();
(EXIT_FAILURE);
no female age specified on the last line of the control file is greater than the last (>0 prob) female age in the age probs file.
maxafage)
<< "ERROR: A female age in " << root + paramfile << " is greater than the maximum allowed in " << root + ageprfile << "." << endl;
<< "--> Largest female age (" << paramfile << ") = " << maxfage << endl;
<< "--> Largest allowable female age (" << ageprfile << ") = " << maxafage << endl;
<< "EXITING...derhrates()..." << endl; cin.get();
(EXIT_FAILURE);

// Derive total # years to process for this arrangement (plus 1). (Includes SMax, but no alpha decision is made at SMax.)
remyrs = max(maxamage-minmage+1, maxafage-minfage+1);
// If the hazard rates file exists use it. Otherwise derive it and write to file (using the details on the last line of the control file).
if (hrexst)
cout << endl << "ALERT: Hazard rates file: " << root + hrfile << " already exists and will be used. (Last line of the control file is ignored.)" << endl;
else
{

42

// Construct the person-specific probability arrays (pmfs and cdfs), shifting each to begin at t=0 which represents the start of retirement.
// Build array of pointers to individual TD probability arrays (start at time t=0 for each person). Then build inner arrays that will
// apply to each individual person retiring, containing remyrs elements. (Do this for both the individual person pmfs and cdfs.)
perspmf = new long double * [npersons];
perscdf = new long double * [npersons];
tdcdf = new long double [remyrs];
hrates = new long double [remyrs];
for (int i=0; i<npersons; ++i)
{
// Check that age values in parameter file are valid. (Invalid gender values are checked in main().)
if (ags[i] < strtage || gndrs[i] == 'M' && ags[i] > maxamage || gndrs[i] == 'F' && ags[i] > maxafage)
{
cout << "ERROR: Invalid age in file control.txt for person #" << i << " of gender=" << gndrs[i] << ": " << ags[i] << endl;
cout << "EXITING...derhrates()..." << endl; cin.get();
exit (EXIT_FAILURE);
}
// Create the individual person PMF and CDF arrays.
perspmf[i] = new long double [remyrs];
perscdf[i] = new long double [remyrs];
// Populate the individual arrays for each person. These are probabilities of death that
// will sum to one and a value exists for each of the remyrs processed by this arrangement.
sumprobs = 0.0;
for (int j=((int) mprobs.size()-1); j>=(ags[i]-strtage); j--)
{
if (gndrs[i] == 'M')
sumprobs = sumprobs + mprobs[j];
else if (gndrs[i] == 'F')
sumprobs = sumprobs + fprobs[j];
}
// Load the individual person PMF array with the applicable conditional PMF value. (Initialize with zeros.)
for (int j=0; j<remyrs; ++j)
perspmf[i][j] = 0.0;
for (int j=0; (gndrs[i] == 'M' && j<(maxamage-(ags[i]-1))) || (gndrs[i] == 'F' && j<(maxafage-(ags[i]-1))); ++j)
{
if (gndrs[i] == 'M')
perspmf[i][j] = mprobs[ags[i]-strtage+j] / sumprobs;
else if (gndrs[i] == 'F')
perspmf[i][j] = fprobs[ags[i]-strtage+j] / sumprobs;
}
// Check that sum of probabilities equals 1 (within floating point precision). If not put out alert.
sumprobs = 0.0;
for (int j=0; j<remyrs; ++j)
sumprobs = sumprobs + perspmf[i][j];
if (sumprobs < (1.00 - pow(0.1,15)) || sumprobs > (1.00 + pow(0.1,15)))
{
cout.setf(ios_base::fixed, ios_base::floatfield); cout.precision(50);
cout << "Alert: Sum of probabilities for Person " << i+1 << " is not 1.00." << endl << "Sum is = " << sumprobs << endl;
}
// Load the individual person cdf array with the applicable cumulative probability.
// And load the TD cdf array with the applicable cumulative probability.
for (int j=0; j<remyrs; ++j)
{
perscdf[i][j] = 0;
for (int k=0; k<=j; ++k)
perscdf[i][j] = perscdf[i][j] + perspmf[i][k];
if (i==0)
tdcdf[j] = perscdf[0][j];

43

else
}

}

tdcdf[j] = tdcdf[j]*perscdf[i][j];

// Build the hazard rates and write to file.
hrates[0] = tdcdf[0];
for (int j=1; j<remyrs; ++j)
hrates[j] = min((tdcdf[j] - tdcdf[j-1]) / (1.0 - tdcdf[j-1]), 1.00);
// Write hazard rate probabilities to file for reference, this file will be read back in if processing years in separate calls.
ofstream hrout(root+hrfile);
hrout.setf(ios_base::fixed, ios_base::floatfield);
for (int j=0; j<remyrs; ++j)
{
hrout.precision(50); hrout << hrates[j] << " "; hrout.precision(0); hrout << "(t=" << j << ")" << endl;
}
hrout.close();

}

}

// Free dynamic memory allocations.
for (int i=0; i<npersons; ++i)
{
delete [] perspmf[i]; perspmf[i] = nullptr;
delete [] perscdf[i]; perscdf[i] = nullptr;
}
delete [] perspmf; perspmf = nullptr;
delete [] perscdf; perscdf = nullptr;
delete [] tdcdf; tdcdf = nullptr;
delete [] hrates; hrates = nullptr;

// Return # years to process for this arrangement, which is 1 less than the # written to the hr file (since SMax is written).
return (remyrs-1);

Code File: getprprobs.cpp
/*
/ Copyright (c) 2014 Chris Rook, All Rights Reserved.
/
/
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
/
/
1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
/
2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the
/
documentation and/or other materials provided with the distribution.
/
3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this
/
software without specific prior written permission.
/
/
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
/
TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
/
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
/
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
/
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
/
EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. (Source: http://opensource.org/licenses/BSD-3-Clause.)
/
/ Filename: getprprobs.cpp
/
/ Function: getprprobs()
/
/ Summary:
/

44

/
This function reads the prior time point's probability file and populates an array with these values. It also returns a vector containing only the
/
bucket #'s from the prior time point that have unique probabilities assigned to them. When deriving the expected probability of ruin for any time
/
point after the next time point it is not necessary to process every bucket if sequential buckets hold the same probability. Therefore the optimize()
/
function only processes those buckets with unique probabilities. Note that the probability of falling into a new larger bucket that was formed from
/
multiple sequential buckets is higher. When at time t=TD-1 or t=SMax-1 (i.e., the first decision point processed via induction) we set V(t,RF(t))=0
/
for all RF(t) > 0. There is no external file to read when processing times t=TD-1 or t=SMax-1.
/
/ Parameters:
/
/
1.) Root directory where input files reside and output files are written.
/
2.) Current year/time point being processed. (These start at 0 and end at one less than the total # being processed by the given arrangement specified
/
in the control file. When we proceed using induction, we start at the end and finish at the beginning.)
/
3.) Total number of years processed by this control file arrangement (i.e., total # of alpha-decisions required).
/
4.) Maximum ruin factor value specified in the control file (i.e., RFMax).
/
5.) Precision value for the ruin factor discretization as specified in the control file (i.e., PR).
/
6.) An empty array of appropriate size (i.e., total # of buckets + 1) long doubles that will be populated by this function.
/
7.) Value of the prior time point's hazard rate (0.0 for fixed TD and t < TD, 1 at TD). This is to determine when the highest possible probability has
/
been reached and all buckets from this time point onward can be treated as one collection.
/
/ Input Files:
/
/
1.) File of probabilities from the prior time point that has already been derived. This file has the form "Year_X_All_Buckets.txt", where X refers to
/
actual time point. Time points start at 0 and end at 1 less than the total # being processed by the given arrangement. Each time point reflects
/
an alpha-decision that needs to be made. Time t=0 reflects the first asset allocation decision at the start of retirement and time T=SMax-1 reflects
/
the last possible asset allocation decision that needs to be made the year before SMax when TD is random. When TD is fixed, TD-1 is the last time
/
point requiring an asset allocation decision. Note that the prior year/time point probabilities being read are from the current value + 1.
/
/ Output Files:
/
/
None.
/
/ Return Value:
/
/
This function returns a vector of bucket #'s with unique probabilities (from the prior timepoint which has already been processed). (It also populates
/
the array Vp whose location is passed to this function.)
/------------------------------------------------------------------------------------------------------------------------------------------------------------------*/
#include "stdafx.h"
vector<long> getprprobs(const string root, const int curyear, const int nyrs, const double rfmax, const int pprec, long double * Vp, const long double prhrate)
{
int inyr, prMax=0;
long * bktPrn, inb;
double inrf, inalpha;
long double inprob, prevprob;
const long nbuckets = (long) (rfmax*pprec + 0.5);
vector<long> PrB;
string ifname;
ifstream fin;
ifname = root + "Year_" + to_string((long long) (curyear+1)) + "_All_Buckets.txt";
// Array to indicate (0/1) where unique probabilities from prior time point reside. By only processing unique probabilities at the prior
// time point we can reduce processing time. (We will always process the first and last buckets.)
bktPrn = new long [nbuckets];
if (curyear == (nyrs-1))
{
for (int i=0; i<nbuckets; ++i)
Vp[i]=0.00;
Vp[nbuckets]=0.00;
}
else {

// If at last decision point (either t=TD-1 or t=SMax-1) then V(t+1,RF(t+1))=0 for all RF(t+1)>0.

// Pruin for RF(t+1) region > RFMax at TD-1/SMax-1.

45

// Open and read the prior time point's probability file and load the values into an array.
fin.open(ifname);
if (fin.is_open())
{
while (!fin.eof())
{
fin >> inyr >> inrf >> inprob >> inalpha;
if (inyr == (curyear + 1))
{
inb = (long) (inrf*pprec + 0.5);
Vp[inb-1]=inprob;
}
}
fin.close();
Vp[nbuckets]=(1.00 - prhrate);
// PRuin for RF(t+1) region > RFMax at all time points except TD-1/SMax-1.
}
else
{
cout << "ERROR: Could not read prior year probability file: " << ifname << endl;
cout << "EXITING...getprprobs()..." << endl; cin.get();
exit (EXIT_FAILURE);
}
}
// Check the data just read in for any issues and load pruning array.
for (long b=1; b <= nbuckets; ++b)
{
// Populate the pruning array with 1's and 0's, where 1 indicates unique probability (last bucket in the sequence).
if (b > 1 && b < nbuckets && Vp[b-1] != Vp[b] && prMax == 0) // Always process first and last buckets individually.
{
if (Vp[b] >= (1.00 - prhrate))
// The next bucket holds the maximum possible probability, stop scanning.
prMax = 1;
bktPrn[b-1] = 1;
// This bucket will account for a unique probability.
}
else if (b > 1 && b < nbuckets)
bktPrn[b-1] = 0;
// This bucket will not account for a unique probability.
else
bktPrn[b-1] = 1;
// First and last buckets will always be processed.
if (b==1)
prevprob=0;

}

// Check for and report issues found in the prior time point's probability file.
if ( (Vp[b-1] < 0) || (Vp[b-1] > (1.00 - prhrate) + (2.0*pow(0.1,16))) || (Vp[b-1] < (prevprob - pow(0.1,15))) )
{
cout.setf(ios_base::fixed, ios_base::floatfield);
cout.precision(50);
cout << "There is an issue with the previous years data (see below):" << endl;
cout << "Vp[" << (b-2) << "] = " << Vp[b-2] << endl;
cout << "Vp[" << (b-1) << "] = " << Vp[b-1] << endl;
cout << "EXITING...getprprobs()..." << endl; cin.get();
exit(EXIT_FAILURE);
}
prevprob = Vp[b-1];

// Build an array with values that are the bucket #'s from the prior time point that have a
// unique probability and thus will be processed during the expected value computation.
for (long b=1; b<=nbuckets; ++b)
if (bktPrn[b-1] == 1)
PrB.push_back(b);

46

// Free memory of dynamically sized arrays.
delete [] bktPrn; bktPrn = nullptr;
// For informational purposes.
cout << "--> # Unique bucket probs at prior time point reported by getprprobs(): " << PrB.size() << endl;

}

// Return the vector PrB.
return PrB;

Code File: optimize.cpp
/*
/ Copyright (c) 2014 Chris Rook, All Rights Reserved.
/
/
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
/
/
1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
/
2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the
/
documentation and/or other materials provided with the distribution.
/
3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this
/
software without specific prior written permission.
/
/
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
/
TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
/
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
/
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
/
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
/
EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. (Source: http://opensource.org/licenses/BSD-3-Clause.)
/
/ Filename: optimize.cpp
/
/ Function: optimize()
/
/ Summary:
/
/
This function derives the optimal alpha and probability values for a single year/time point and bucket range specified by the call. When the first
/
and last buckets are set to 0 this function makes a single pass over the given year/time point and approximates the bucket number where heavy algorithm
/
pruning begins. Heavy pruning sets alpha=pa which is a decision rule that takes effect when the optimal probability of ruin exceeds the first k
/
decimal places of the maximum possible probability of ruin for this year/arrangment, where k is the last parameter specified on the first line of the
/
control file. Thus a value of 4.00 indicates that pruning begins when the probability of ruin exceeds the first 4 decimal places of the highest possible
/
probability of ruin for this year/arrangement. Note that this value is stored locally in the variable prnpwr and only has an impact when TD is random.
/
The initial pass approximates this value by only considering alpha=pa. Optimal bucket sizes are then formed based off of this approximation and then
/
the actual calls to this function are made where the point of heavy pruning is determined exactly. The approximation is only used to estimate optimal
/
bucket collection sizes for threaded processing. For all other cases the specific buckets listed are processed and written to a file with the bucket
/
boundaries contained in the file name. The function combfiles() is invoked by main() after all buckets for a year/time point have been processed and
/
written to a file and it aggregates all individual files then deletes them leaving a single file containing all relevant optimal values for each year/
/
time point. Once the final year/time point has been processed (i.e., t=0) the individual files for each year/time point are aggregated into final
/
(horizontal & vertical) files for the specific arrangement being modeled and the year-specific files are deleted.
/
/ Parameters:
/
/
1.) Root directory where input files reside and output files are written.
/
2.) Year value to process. (This is the timepoint value which always starts at zero and ends at # time points to process for this arrangement less 1.)
/
3.) Array containing all values from the first line of the control file.
/
4.) Array containing the precision values from the 2nd line of the control file.
/
5.) Array containing the bucket limits for this call.
/
6.) Array containing the hazard rates derived for this arrangement.
/
7.) Array of optimal probabilities derived for the prior year/time point (i.e., current year/time point + 1)
/
8.) Vector of unique bucket endpoints for prior year where collections of buckets with equal probabilities are treated as a single bucket in the
/
expected value calculation.

47

/
/ Input Files:
/
/
None. (The file of prior probabilities is read by the function getprprobs() and passed to this function via the parameter Vp.)
/
/ Output Files:
/
/
1.) The file "Year_Y_Buckets_S_thru_E.txt" is written to directory root. (Where Y is the time point, S is the start bucket #, and E is the end bucket #.
/
Once processing for this year has ended these files are aggregated by combfiles() into a single file named "Year_Y_All_Buckets.txt".)
/
/ Return Value:
/
/
The approximate bucket # where heavy algorithm pruning begins when invoking this function with start and end buckets equal to 0. Otherwise return
/
the value 0.
/-------------------------------------------------------------------------------------------------------------------------------------------------------------*/
#include "stdafx.h"
int optimize(const string root, const int y, const double prms[8], const int prc[2], const long bkts[2], const long double * sprobs, const long double * Vp, const
vector<long> & PrB)
{
const int pr=prc[0], pa=prc[1];
const double smn=prms[0], svr=prms[1], bmn=prms[2], bvr=prms[3], cv=prms[4], rfmax=prms[5], er=prms[6], prnpwr=prms[7];
const long nbuckets=(long) (rfmax*pr + 0.5), nuqbkts=(long) PrB.size();
int prunepnt = 0, ties, stalpha=0;
long stbkt=bkts[0], endbkt=bkts[1];
double rf, alpha, OPT_alpha, mean, std, * Ac;
long double OPT_pruin = 1.0, cdfval, eprob, * Vc, rhs_cdf, lhs_cdf, pruin, tiethresh, pruneprob;
boost::math::normal normdist;
ofstream fout;
// Create the arrays to hold current year optimal alphas and probabilities.
Ac = new double [nbuckets];
Vc = new long double [nbuckets];
// Before processing the current year, run through all buckets using only alpha=100% to find the pruning point.
if (stbkt == 0 && endbkt == 0)
{
stalpha = pa;
stbkt = 1;
endbkt = nbuckets;
}
// One-half of the maximum possible PRuin is the tie threshold, in general.
tiethresh = (0.5)*((1.00 - sprobs[y]));
// Prune point #2 begins at set number of decimals after max probability for this arrangement.
pruneprob = floor(pow(10.0,prnpwr)*(1.00 - sprobs[y])) / pow(10.0,prnpwr);
if (bkts[0] == 0 && bkts[1] == 0)
{
cout.setf(ios_base::fixed, ios_base::floatfield);
cout.precision(10); cout << "--> Value of pruneprob reported by optimize() is: " << pruneprob << endl;
}
for (long b = stbkt; b <= endbkt; ++b)
// PRuin at any future time point.
{
ties = 0;
// Initialize for le/gt comparisons.
rf = (double) b / pr;
// Derive ruin factor.
OPT_alpha = 99.00;
// Initialize to unrealistic value.
OPT_pruin = 99.00;
// Initialize to unrealistic value.
for (int a=stalpha; a <= pa; ++a)
{

48

if ((prunepnt == 0 && (a == stalpha || OPT_pruin > 0.00)) || (prunepnt == 1 && a == pa))
{
alpha = (double) a / pa;
mean = (1.00-er)*(1.00 + alpha*smn + (1.00-alpha)*bmn);
std = (1.00-er)*sqrt(pow(alpha,2)*svr + pow(1.00-alpha,2)*bvr + 2.00*(alpha)*(1.00-alpha)*(cv));
normdist = boost::math::normal(mean,std);
cdfval = cdf(normdist, rf);
if (cdfval == 1.00)
eprob = Vp[nbuckets];
else
{
//----------------------------------------------------------------------------------------------------------------------------//
rhs_cdf = 1.00;
lhs_cdf = cdf(normdist, rf*(1 + (pr/1.5)));
eprob = (rhs_cdf - lhs_cdf)*Vp[0];
// First bucket, unique processing.
rhs_cdf = lhs_cdf;
for (long pb=2; pb <= nuqbkts; ++pb)
// All others but last bucket, std proc for unq probs only.
{
lhs_cdf = cdf(normdist, rf*(1.0 + pr/(PrB[pb-1] + 0.5)));
eprob = eprob + (rhs_cdf - lhs_cdf)*Vp[PrB[pb-1]-1];
rhs_cdf = lhs_cdf;
}
eprob = eprob + (rhs_cdf - cdfval)*Vp[nbuckets];
// Last bucket, unique processing.
eprob = eprob / (1.00 - cdfval);
// Make the probability conditional.
//----------------------------------------------------------------------------------------------------------------------------//
}
// Deal with numerical instability near zero.
if (ties == 0)
{
pruin = (1.00 - sprobs[y])*(cdfval + eprob - (cdfval*eprob));
if (pruin > tiethresh)
ties = 1;
}
// Deal with numerical instability near one.
if (ties == 1)
pruin = (1.00) - (sprobs[y] + (1.00-cdfval)*(1.00-eprob) - sprobs[y]*(1.00-cdfval)*(1.00-eprob));
// Update optimal values.
if ((a == 0) || (ties == 0 && pruin < OPT_pruin) || (ties == 1 && pruin <= OPT_pruin))
{
OPT_pruin = pruin;
OPT_alpha = alpha;
}

}
}
// Load optimal values into correspnding arrays.
Ac[b-1] = OPT_alpha;
Vc[b-1] = OPT_pruin;

// Set pruning for this timepoints processing, allow for floating point precision diffs.
if (prunepnt == 0 && Vc[b-1] >= (pruneprob - (pow(0.1,16) + pow(0.1,17))))
{
prunepnt = 1;
// Initial run finds the (approximate) bucket where pruning starts then uses this information to split the remaining buckets by
// sizes that will run in shortest time. For this run, end the call as soon as the bucket number has been found.
if (bkts[0] == 0 && bkts[1] == 0)
{
cout << "--> Pruning RC returned by optimize() is: " << b << endl;

49

}

}

}

return b;

// Write current year to file with bucket boundaries contained in the name.
string outfile = root + "Year_" + to_string((long long) y) + "_Buckets_" + to_string((long long) stbkt) + "_thru_" + to_string((long long) endbkt) + ".txt";
fout.open(outfile);
for (long b = stbkt; b <= endbkt; ++b)
{
fout.setf(ios_base::fixed, ios_base::floatfield);
fout << y << " ";
fout.precision(10); fout << (double) b / pr << " ";
fout.precision(50); fout << Vc[b-1] << " ";
fout.precision(10); fout << Ac[b-1] << endl;
}
fout.close();

}

// Free dynamic memory allocations and exit.
delete [] Ac; Ac = nullptr;
delete [] Vc; Vc = nullptr;
return 0;

Code File: combfiles.cpp
/*
/ Copyright (c) 2014 Chris Rook, All Rights Reserved.
/
/
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
/
/
1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
/
2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the
/
documentation and/or other materials provided with the distribution.
/
3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this
/
software without specific prior written permission.
/
/
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
/
TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
/
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
/
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
/
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
/
EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. (Source: http://opensource.org/licenses/BSD-3-Clause.)
/
/ Filename: combfiles.cpp
/
/ Function: combfiles()
/
/ Summary:
/
/
Function to aggregate the individual bucket-specific output files written by optimize() into single files for each time point. These files are
/
subsequently read by getprprobs() when processing the next time point. At the last time point (i.e., year=0) all files for the individual time
/
points are aggregated into 3 files then deleted. The first is a (space-delimited) vertical file containing all data for all time points. This
/
file is named FinalResults_V.txt. The remaining 2 files are comma-separated horizontal files for import into a spreadsheet, and these files are
/
named FinalAlphaResults_H.csv and FinalProbResults_H.csv.
/
/ Parameters:
/
/
1.) Root directory where input files reside and output files are written.
/
2.) Current year/time point number being processed. (Years start at 0 and end at one less than the total # of years being processed by the given

50

/
arrangement specified in the control file.)
/
3.) Number of concurrent processes used to derive the optimal values for each time point.
/
4.) Number of buckets processed per thread within each time point. (Except the last one. This is specific to each time point.)
/
5.) Maximum ruin factor value specified in the control file (i.e., RFMax).
/
6.) Precision value for the ruin factor discretization as specified in the control file (i.e., PR).
/
7.) Total number of years processed by the arrangement specified in the control file.
/
/ Input Files:
/
/
1.) For a given time point, input files of the form "Year_X_Buckets_S_thru_E.txt" are read from the root directory. (Where X is the time point,
/
S is the start bucket #, and E is the end bucket #. This function aggregates these individual files into a single file named: Year_X_All_Buckets.txt.
/
2.) Year_X_All_buckets.txt (Written by this function for each time point and all such files are read by this function and aggregated at the last time point,
/
which is year=0.)
/
/ Output Files:
/
/
1.) Year_X_All_Buckets.txt, where X reflects the year/time point being processed.
/
2.) FinalResults_V.txt (Vertical text file of all results from all years processed.)
/
3.) FinalAlphaResults_H.csv (Horizontal csv file of just optimal alphas. Use this file to simulate and confirm the result.)
/
4.) FinalProbResults_H.csv (Horizontal csv file of optimal probabilities.)
/
/
Note: The last 3 output files above are only written at the last time point when all processing ends (i.e., year=0). The first file is written at every
/
time point.
/
/ Return Value:
/
/
None.
/---------------------------------------------------------------------------------------------------------------------------------------------------------------------*/
#include "stdafx.h"
void combfiles(const string root, const int curyear, const int pllruns, const long bcktsprun, const double rfmax, const int pprec, const int nyrs)
{
int numfiles=0, fnlyear;
const long nbkts=(long) (rfmax*pprec + 0.5);
long fnlbkt;
double fnlrf, fnlalpha;
long double fnlprob;
string fndfname, fndfile, concatyr, concatall;
ifstream chk4file, normin;
ofstream fnloutp, fnlouta;
numfiles = 0;
for (int j=1; j<=pllruns; ++j)
{
if (j <= (pllruns-1))
fndfname = "Year_" + to_string((long long) curyear) + "_Buckets_" + to_string((long long) bcktsprun*(j-1) + 1) + "_thru_"
+ to_string((long long) bcktsprun*(j)) + ".txt";
else
fndfname = "Year_" + to_string((long long) (curyear)) + "_Buckets_" + to_string((long long) bcktsprun*(j-1) + 1) + "_thru_"
+ to_string((long long) nbkts) + ".txt";
fndfile = root + fndfname;
// Construct string for concatenation.
concatyr=concatyr + " " + fndfile;

}

// Attempt to open file and update file counter if successful.
chk4file.open(fndfile);
if (chk4file.is_open())
numfiles = numfiles + 1;
chk4file.close();

51

if (numfiles == pllruns)
{
// Concatenate files into 1 and remove individual files.
concatyr="type" + concatyr + " > " + root + "Year_" + to_string((long long) curyear) + "_All_Buckets.txt 2>nul & del " + root + "Year_"
+ to_string((long long) (curyear)) + "_Buckets*.txt";
system(concatyr.c_str());
}
else
{
cout << "ERROR: Attempt to combine files after year=" << curyear << " has failed." << endl;
cout << "EXITING...combfiles()..." << endl; cin.get();
exit (EXIT_FAILURE);
}
// If at last time point then build the normalized and transposed final data files and delete all intermediate data files.
if (curyear == 0)
{
// Build normalized data file.
concatall="type " + root + "Year_" + to_string((long long) curyear) + "_All_Buckets.txt";
for (int j=curyear+1; j<nyrs; ++j)
concatall = concatall + " " + root + "Year_" + to_string((long long) j) + "_All_Buckets.txt";
concatall=concatall + " > " + root + "FinalResults_V.txt 2>nul & del " + root + "Year_*.txt";
system(concatall.c_str());
// Build transposed probability/alpha data files.
long double * * fnlprobs;
double * * fnlalphas;
fnlprobs = new long double * [nyrs];
fnlalphas = new double * [nyrs];
for (int j=0; j<nyrs; ++j)
{
fnlprobs[j] = new long double [nbkts];
fnlalphas[j] = new double [nbkts];
}

//
//
//
//
//

Pointer to long double pointer for probabilities.
Pointer to double pointer for alphas.
Use long double pointer to hold array of long double pointers.
Use double pointer to hold array of double pointers.
(A pointer to any type can hold an array of that type.)

// Each array value now points to array of long doubles (2D array).
// Each array value now points to array of doubles (2D array).

// Open normalized file just created, transpose it and write to csv file.
normin.open(root+"FinalResults_V.txt");
if (normin.is_open())
{
while (!normin.eof())
{
normin >> fnlyear >> fnlrf >> fnlprob >> fnlalpha;
fnlbkt = (long) (fnlrf*pprec + 0.5);
if (fnlbkt >= 1 && fnlbkt <= (long) (rfmax*pprec + 0.5))
{
fnlprobs[fnlyear][fnlbkt-1] = fnlprob;
fnlalphas[fnlyear][fnlbkt-1] = fnlalpha;
}
}
normin.close();
}
else
{
cout << "ERROR: Could not read file: " << root+"FinalResults_V.txt" << endl;
cout << "EXITING...combfiles()..." << endl; cin.get();
exit (EXIT_FAILURE);
}
// Transpose and write to file.
fnloutp.open(root+"FinalProbResults_H.csv"); fnloutp.setf(ios_base::fixed, ios_base::floatfield); fnloutp << "RF";
fnlouta.open(root+"FinalAlphaResults_H.csv"); fnlouta.setf(ios_base::fixed, ios_base::floatfield); fnlouta << "RF";

52

for (long i=0; i<=nbkts; ++i)
// The first line (i=0) of the csv files holds column headers.
{
if (i > 0)
{
float _rf_; _rf_ = (float) i/pprec;
fnloutp.precision(10); fnloutp << _rf_;
fnlouta.precision(10); fnlouta << _rf_;
}
for (int j=0; j<nyrs; ++j)
{
if (i == 0)
{
fnloutp << ", Time (t=" << j << ")"; fnlouta << ", Time (t=" << j << ")";
}
else
{
fnloutp.precision(50); fnloutp << "," << fnlprobs[j][i-1];
fnlouta << "," << fnlalphas[j][i-1];
}
}
fnloutp << endl; fnlouta << endl;
}
fnloutp.close();
fnlouta.close();

}

}

// Free dynamic memory allocations.
for (int j=0; j<(nyrs); ++j)
{
delete [] fnlprobs[j]; fnlprobs[j] = nullptr;
delete [] fnlalphas[j]; fnlalphas[j] = nullptr;
}
delete [] fnlprobs; fnlprobs = nullptr;
delete [] fnlalphas; fnlalphas = nullptr;

Code File: stdafx.h
/*
/ Copyright (c) 2014 Chris Rook, All Rights Reserved.
/
/ Filename: stdafx.h (header file to include)
/ ----------------------------------------------------*/
#pragma once
// Header files
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include

"targetver.h"
<stdio.h>
<tchar.h>
<iostream>
<fstream>
<string>
<math.h>
<stdlib.h>
<boost/math/distributions.hpp>
<boost/algorithm/string.hpp>
<boost/thread/thread.hpp>
<windows.h>
<dos.h>
<vector>

53

using namespace std;
// Define constants.
const string paramfile="control.txt";
const string ageprfile="ageprobs.txt";
const string hrfile="hrates.txt";

// Name of parameter control input file.
// Name of age probability input file.
// Name of hazard rate probability output file.

// Define Function prototypes.
int optimize(const string root, const int strtyr, const double prms[8], const int prc[2], const long bkts[2], const long double * sprobs, const long double * Vp, const
vector<long> & PrB);
int derhrates(const string root, const int numpersons, const char * ppl, const int * ags, const bool hrexst);
void combfiles(const string root, const int curyear, const int pllruns, const long bcktsprun, const double rfmax, const int pprec, const int nyears);
vector<long> getprprobs(const string root, const int curyear, const int nyrs, const double rfmax, const int pprec, long double * Vp, const long double prhrate);

54

Output File: hrates.txt
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

0.00000000000000031197163361710200000000000000000000
0.00000000000021922536698186529000000000000000000000
0.00000000001136915190873025700000000000000000000000
0.00000000020003863625408732000000000000000000000000
0.00000000195081090367563750000000000000000000000000
0.00000001312223201754271500000000000000000000000000
0.00000006827600465499560600000000000000000000000000
0.00000029370323648036120000000000000000000000000000
0.00000109220107904646070000000000000000000000000000
0.00000361453274149872190000000000000000000000000000
0.00001084016753252347100000000000000000000000000000
0.00002984584354025825700000000000000000000000000000
0.00007620909696035615500000000000000000000000000000
0.00018188647334656542000000000000000000000000000000
0.00040760821017616586000000000000000000000000000000
0.00085904264210663790000000000000000000000000000000
0.00170489014561057340000000000000000000000000000000
0.00319359439131450570000000000000000000000000000000
0.00565897467490714380000000000000000000000000000000
0.00949932496688625510000000000000000000000000000000
0.01511818909158774700000000000000000000000000000000
0.02287508163829699400000000000000000000000000000000
0.03301592282445740800000000000000000000000000000000
0.04563765178543890300000000000000000000000000000000
0.06069174610916736000000000000000000000000000000000
0.07796388736906334500000000000000000000000000000000
0.09713739625524489600000000000000000000000000000000
0.11784066318685273000000000000000000000000000000000
0.13957595061039985000000000000000000000000000000000
0.16190010225257540000000000000000000000000000000000
0.18446735111765136000000000000000000000000000000000
0.20713558135530538000000000000000000000000000000000
0.22917196829287370000000000000000000000000000000000
0.25057440720421520000000000000000000000000000000000
0.27097719546324184000000000000000000000000000000000
0.29008553190192188000000000000000000000000000000000
0.30800961581846242000000000000000000000000000000000
0.32507762001542762000000000000000000000000000000000
0.34246475410896637000000000000000000000000000000000
0.36010444680717479000000000000000000000000000000000
0.37955393469636328000000000000000000000000000000000
0.39930708184585023000000000000000000000000000000000
0.42023780462667704000000000000000000000000000000000
0.44104841116863625000000000000000000000000000000000
0.46608832292708718000000000000000000000000000000000
0.49521022619328664000000000000000000000000000000000
0.53249502096344137000000000000000000000000000000000
0.54626727760656602000000000000000000000000000000000
0.59866287346598623000000000000000000000000000000000
0.75414991674801224000000000000000000000000000000000
0.99999999997023048000000000000000000000000000000000

(t=0)
(t=1)
(t=2)
(t=3)
(t=4)
(t=5)
(t=6)
(t=7)
(t=8)
(t=9)
(t=10)
(t=11)
(t=12)
(t=13)
(t=14)
(t=15)
(t=16)
(t=17)
(t=18)
(t=19)
(t=20)
(t=21)
(t=22)
(t=23)
(t=24)
(t=25)
(t=26)
(t=27)
(t=28)
(t=29)
(t=30)
(t=31)
(t=32)
(t=33)
(t=34)
(t=35)
(t=36)
(t=37)
(t=38)
(t=39)
(t=40)
(t=41)
(t=42)
(t=43)
(t=44)
(t=45)
(t=46)
(t=47)
(t=48)
(t=49)
(t=50)

Note: “hrates.txt” is a file derived by the function derhrates() and is written for random TD
arrangements. It contains the discrete hazard rates for the person/MPU defined on the last line
of the control file. (See Section 4.6.2 for reference.) The probabilities shown above are from the
5th “control.txt” example (MPU of size N=9). All random arrangements have a single set of
discrete hazard rates which is the only random-specific input required by the random TD model.
Because an MPU is defined by a single hazard rate file the MPU size does not impact runtime.
In the example above, SMax=50 and a maximum of 50 decisions need to be made, the first at t=0
and the last at t=SMax-1=49. No decision is made at time t=SMax since P(TD > SMax) = 0.
55

Output File: FinalResults_V.txt
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

0
0
0
0

0.0016000000
0.0017000000
0.0018000000
0.0019000000

0.00000000000000000001927209692083706200000000000000
0.00000000000000000005785455286366354000000000000000
0.00000000000000000016213529004089414000000000000000
0.00000000000000000042735423336906475000000000000000

0.1980000000
0.1990000000
0.1990000000
0.2000000000

.
.
.

10
10
10
10

0.0108000000
0.0109000000
0.0110000000
0.0111000000

0.00000000222864883182142450000000000000000000000000
0.00000000256035079632214650000000000000000000000000
0.00000000293658568468094670000000000000000000000000
0.00000000336267466389999980000000000000000000000000

0.2140000000
0.2140000000
0.2150000000
0.2150000000

.
.
.

20
20
20
20

0.0296000000
0.0297000000
0.0298000000
0.0299000000

0.00001079481991831277300000000000000000000000000000
0.00001125172852740063400000000000000000000000000000
0.00001172534597344354000000000000000000000000000000
0.00001221616101523663900000000000000000000000000000

0.2360000000
0.2360000000
0.2370000000
0.2370000000

.
.
.

30
30
30
30

0.0934000000
0.0935000000
0.0936000000
0.0937000000

0.00775593991107923190000000000000000000000000000000
0.00780733692967890570000000000000000000000000000000
0.00785894808515004080000000000000000000000000000000
0.00791077179182503140000000000000000000000000000000

0.3640000000
0.3640000000
0.3650000000
0.3650000000

.
.
.

40
40
40
40

0.0001000000
0.0002000000
0.0003000000
0.0004000000

0.00000000000000000000000000000000000000085482880000
0.00000000000000000000000000000000000000088469110000
0.00000000000000000000000000000000000000091611610000
0.00000000000000000000000000000000000000094839560000

0.1440000000
0.1440000000
0.1440000000
0.1440000000

.
.
.

45
45
45
45

0.4978000000
0.4979000000
0.4980000000
0.4981000000

0.08876385008462925300000000000000000000000000000000
0.08883038694792594300000000000000000000000000000000
0.08889690950186190900000000000000000000000000000000
0.08896342425482126100000000000000000000000000000000

0.4170000000
0.4180000000
0.4200000000
0.4210000000

.
.
.

46
46
46
46

0.5223000000
0.5224000000
0.5225000000
0.5226000000

0.04775335689522756900000000000000000000000000000000
0.04780431578767197600000000000000000000000000000000
0.04785521417275828700000000000000000000000000000000
0.04790604639994011600000000000000000000000000000000

0.8510000000
0.8530000000
0.8550000000
0.8560000000

.
.
.

47
47
47
47

0.8091000000
0.8092000000
0.8093000000
0.8094000000

0.00104883991561901090000000000000000000000000000000
0.00105293787325091270000000000000000000000000000000
0.00105704928401912120000000000000000000000000000000
0.00106117491595328910000000000000000000000000000000

0.1740000000
0.1740000000
0.1750000000
0.1750000000

Note: Select rows from the file “FinalResults_V.txt” are displayed. This is a vertical file written
by the application containing all model output. The values shown are from the 6th “control.txt”
example. Column definitions are: time point [t], ruin factor RF(t), minimum probability of ruin
VR(t, RF(t)), and optimal asset allocation αR(t, RF(t)). In addition, there are 2 horizontal files
which together contain the same information as the single vertical file and these are named:
“FinalAlphaResults_H.csv” & “FinalProbResults_H.csv”. To simulate and confirm the result
use the file “FinalAlphaResults_H.csv”. A method for doing this is discussed in the next section.
56

Additional Notes:


The source code provided has not been independently validated. The user should conduct
full code validation before using it within a production environment.7



To simulate a specific solution load the file “FinalAlphaResults_H.csv” into Microsoft®
Excel® and use a simulator such as Oracle® Crystal Ball®. A vlookup table is used to
retrieve the optimal alpha at each time point.8 (Recall that each RF(t) centers its bucket.)



The basic approach for coding the DP suggested in Figure 5 was altered to avoid
unacceptable runtimes. Most coding effort was spent modifying the implementation to
reduce runtimes which also has the undesirable effect of complicating the code.



This implementation has been tested on the Windows 7 Professional and Windows 7
Home Edition operating systems. Some file management operations in the code are
Windows-specific and would need to be updated for use on other operating systems.



The code provided will produce a similar but slightly better solution than was presented
in Figures 6, 7, and 8. Runtimes are now such that there is no longer a need to partition
the RF(t) dimension when discretizing it as was done in those solutions.



Increasing the RF(t) discretization precision PR results in a more accurate numerical
approximation whereas increasing the alpha precision Pα lowers the probability of ruin
V(t, RF(t)).

***

7
8

The author would greatly appreciate being informed of code errors/bugs as well as techniques to reduce runtimes.
A template for simulating the optimal solution will be supplied upon request.

57

REFERENCES
Anton, Howard, 1988, Calculus with Analytic Geometry 3rd Edition (John Wiley & Sons, New
York, NY).
Bengen, William P., 1994, Determining withdrawal rates using historical data, Journal of
Financial Planning 7(4), 171-180.
Bengen, William P., 2006, Baking a withdrawal plan 'layer cake' for your retirement clients,
Journal of Financial Planning 19(8), 44-51.
Blanchett, David M., 2007, Dynamic allocation strategies for distribution portfolios:
Determining the optimal distribution glide-path, Journal of Financial Planning 20(12), 68-81.
Box, George E.P., Gwilym M. Jenkins, and Gregory C. Reinsel, 1994, Time Series Analysis
Forecasting and Control 3rd Edition (Prentice Hall, Englewood Cliffs, NJ).
Bradley, Stephen P., Arnoldo C. Hax, and Thomas L. Magnanti, 1977, Applied Mathematical
Programming (Addison-Wesley Publishing Company, Reading, MA).
Brigham, Eugene F., Michael C. Ehrhardt, 2008, Financial Management Theory and Practice
12th Edition (South-Western CENGAGE Learning, Mason, OH).
Charlson, Josh, Michael Herbst, Kailin Liu, Laura P. Lutton, and John Rekenthaler, 2009,
Target-date series research paper: Industry survey (Morningstar, Inc., Chicago, IL).
Cohen, J., Grant Gardner, and Yuan-An Fan, 2010, The date debate: Should target-date fund
glide-paths be managed "to" or "through" retirement? (Russell Research, Russell Investments,
Seattle, WA).
Cooley, Philip L., Carl M. Hubbard, and Daniel T. Walz, 1998, Retirement savings: Choosing a
withdrawal rate that is sustainable, American Association of Individual Investors Journal 10(3),
16–21.
Cooley, Philip L., Carl M. Hubbard, and Daniel T. Walz, 2011, Portfolio success rates: Where to
draw the line, Journal of Financial Planning 24(4), 48-60.
Damodaran, Aswath, Financial Database of Historical Stock and Bond Returns (New York
University, New York) <http://pages.stern.nyu.edu/~adamodar/>
Fan, Yuan A., Steve Murray and Sam Pittman, 2013, Optimizing Retirement Income: An
Adaptive Approach Based on Assets and Liabilities, The Journal of Retirement 1(1), 124-135.
The Federal Reserve Bank of Minneapolis, CPI-U Data 1913 – Present,
<http://www.minneapolisfed.org/community_education/teacher/calc/hist1913.cfm>

58

Finke, Michael, Wade D. Pfau, and David M. Blanchett, 2013, The 4 percent rule is not safe in a
low-yield world, Journal of Financial Planning 26(6), 46-55.
Guyton, Jonathan T., 2004, Decision rules and portfolio management for retirees: Is the 'safe'
initial withdrawal rate too safe? Journal of Financial Planning 17(10), 54-62.
ING Investment Management and ING Retirement Research Institute, 2012, Participant
preferences in target-date funds: Examining perceptions and expectations among target-date
users and non-users, <http://graphics8.nytimes.com/packages/pdf/tdf_white_paper.pdf>
Irlam, Gordon, and Joseph Tomlinson, 2014, Retirement income research: What can we learn
from economics? Journal of Retirement 1(4), 118-128.
Kitces, Michael E., 2012, Utility Functions in Financial Planning - A New Framework For
Decision Making? Nerd's Eye View, April 3, <http://www.kitces.com/blog/utility-functions-infinancial-planning-a-new-framework-for-decision-making/>
Kotlikoff, Laurence J., 2008, Economics' approach to financial planning. Journal of Financial
Planning 21(3), 42-52.
Markowitz, Harry, 1952, Portfolio selection, The Journal of Finance 7(1), 77–91.
Milevsky, Moshe A., and Huaxiong Huang, 2010, Spending retirement on planet vulcan: The
impact of longevity risk aversion on optimal withdrawal rates, Financial Analysts Journal 67(2),
45-58.
Milevsky, Moshe A., and Chris Robinson, 2005, A sustainable spending rate without simulation,
Financial Analysts Journal 61(6), 89-100.
Mitchell, Olivia S., and Kent Smetters, 2013, The Market for Retirement Financial Advice
(Pension Research Council) (Oxford University Press, Oxford).
Pfau, Wade D., and Michael E. Kitces, 2014, Reducing retirement risk with a rising equity glidepath. Journal of Financial Planning 27(1), 38-45.
Pye, Gordon B., 2000, Sustainable investment withdrawals, Journal of Portfolio Management
26(4), 73-83.
Ross, Sheldon M., 2007, Introduction to Probability Models 9th Edition (Academic Press, New
York, NY).
Scott, Jason S., William F. Sharpe, and John G. Watson, The 4% rule - At what price? Journal of
Investment Management 7(3), 31-48.
Stout, R. Gene, 2008, Stochastic optimization of retirement portfolio asset allocations and
withdrawals. Financial Services Review 17(1), 1-15.
59

