Modeling A Systematic Withdrawal Plan: A Stochastic Algorithm to Estimate Initial Fund Requirement

Abstract

A systematic withdrawal plan refers to monthly withdrawals from an appreciating fund, especially to provide retirement income. This study proposes a stochastic algorithm that outputs the initial amount required to sustain a systematic withdrawal plan for a target duration with a 90% safety margin, given an initial monthly withdrawal amount. The systematic withdrawal plan is based on the S&P 500. Initially, a deterministic function modeling a systematic withdrawal plan was created. This study showed that the monthly growth rate of the S&P 500 is normally distributed with a mean of 0.6933% per month and a standard deviation of 4.0266% per month. By randomly sampling monthly growth rate values from the distribution, stochastic simulations were run to identify the distribution of the fund’s longevity. The Kolmogorov–Smirnov test suggested that the simulated values of the fund’s longevity were log-normally distributed. Initially, the algorithm uses the deterministic function to find the 50% safety margin value for the initial fund requirement. Then, the algorithm uses stochastic simulations to find the log-normal longevity distribution. The 90th percentile value of the log-normal longevity distribution is substituted into the deterministic function to calculate the 90% safety margin value for the initial fund requirement of a systematic withdrawal plan.

Share and Cite:

Dhamane, S. (2025) Modeling A Systematic Withdrawal Plan: A Stochastic Algorithm to Estimate Initial Fund Requirement. Journal of Mathematical Finance, 15, 155-168. doi: 10.4236/jmf.2025.151007.

1. Introduction

To oppose inflation’s depreciating effect on retirement savings, individuals seek to invest their savings to ensure appreciation over time. A systematic withdrawal plan refers to monthly withdrawals from an appreciating fund, especially to provide retirement income [1]. After a certain duration, the increasing withdrawals due to inflation will surpass the income earned by the appreciation rate of the retirement fund. This will cause the fund’s value to reach zero. Hence, the aim of this study is to develop an algorithm to find the initial amount required to sustain a systematic withdrawal plan invested in the S&P 500 for a target duration, given an initial monthly withdrawal amount. Since the S&P 500 is widely diversified across several industries, individuals seek to invest their retirement funds in S&P 500 index funds to minimize risk caused by random chance [2]. Nevertheless, the S&P 500 displays randomness. To model systematic withdrawal plans accurately, it is vital to investigate if this randomness in the monthly growth rate of the S&P 500, as shown in Figure 1, is stochastic. Moreover, systematic withdrawal plans are more vulnerable to randomness in the S&P 500 than regular investments, as low growth rates in the initial months can cause the fund to run out quicker than expected since money is withdrawn from the fund periodically. Hence, to account for the higher risk posed by systematic withdrawal plans, a large safety margin must be determined by comparing a deterministic function modeling the fund with the stochastic simulations modeling the fund. Inflation rate fluctuations in the USA are trivial compared to the S&P 500 fluctuations; hence, the monthly inflation rate is assumed to be constant at 1.0021 (Average monthly inflation rate from 2010 to 2024) [3]. Section 2 involves deriving a deterministic function to model the systematic withdrawal plan, which uses the mean monthly growth rate of the S&P 500 from 2010 to 2024 as a constant. This function is used in the stochastic algorithm in Section 3. Section 3 involves exploratory data analysis to create a log-normal distribution of the fund’s longevity. The stochastic algorithm in this study is a handy tool for estimating the initial fund requirement with 90% certainty based on the target duration and withdrawal amounts.

Figure 1. Monthly growth rate of the S&P 500 (2010-2024).

2. Deterministic Function to Model a Systematic Withdrawal Plan

2.1. Variables

We derive a deterministic function to model a systematic withdrawal plan. Table 1 outlines the variables used in this function.

Table 1. Definitions of the variables.

Variable

Definition

α

The initial amount in the fund

rg

The monthly growth rate of the S&P 500

β

The initial monthly withdrawal amount

ri

The monthly inflation rate in the United States

2.2. Deriving Expressions for Months 1 and 2

Procedure for Month 1:

Step 1: We multiply the initial amount in the fund α by ( 1+ r g / 100 ) to calculate the appreciated amount at Month 1’s end:

α( 1+ r g 100 )

Step 2: We multiply the initial monthly withdrawal amount β by ( 1+ r i / 100 ) to calculate the amount to be withdrawn at Month 1’s end:

β( 1+ r i 100 )

Step 3: We subtract the amount to be withdrawn at Month 1’s end from the appreciated amount at Month 1’s end:

f( 1 )=α( 1+ r g 100 )β( 1+ r i 100 ) (1)

Step 4: We set the amount at Month 1’s end to the amount at the start of Month 2.

Procedure for Month 2:

Step 1: We multiply the amount at the start of Month 2 by ( 1+ r g / 100 ) to calculate the appreciated amount at Month 2’s end:

( α( 1+ r g 100 )β( 1+ r i 100 ) )( 1+ r g 100 )

Step 2: We multiply the amount to be withdrawn at Month 1’s end by ( 1+ r i / 100 ) to calculate the amount to be withdrawn at Month 2’s end:

β ( 1+ r i 100 ) 2

Step 3: We subtract the amount to be withdrawn at Month 2’s end from the appreciated amount at Month 2’s end, which gives

f( 2 )=( α( 1+ r g 100 )β( 1+ r i 100 ) )( 1+ r g 100 )β ( 1+ r i 100 ) 2 =α ( 1+ r g 100 ) 2 ( 1+ r g 100 )β( 1+ r i 100 )β ( 1+ r i 100 ) 2 (2)

Step 4: We set the amount at Month 2’s end to the amount at the start of Month 3. At this stage, two simplifications avoid excessive complexity in the model equation.

Simplification 1: A variable mg is set to the value ( 1+ r g 100 ) .

Simplification 2: A variable mi is set to the value ( 1+ r i 100 ) .

Using the simplifications, the expression for the fund’s value at Month 1’s end is:

f( 1 )=α m g β m i (3)

Using the simplifications, the expression for the fund’s value at Month 2’s end is:

f( 2 )=α m g 2 m g β m i β m i 2 (4)

The procedure used for Months 1 and 2 is repeated for Months 3, 4, 5, and so on until the fund reaches zero.

2.3. Deriving Expressions for Every Month

Table 2 illustrates the expression for the fund’s value at every month’s end after withdrawal. The first five months have been included in the table for representation purposes.

Table 2. The amount in the fund at the month’s end after withdrawal.

Month

The amount in the fund at the month’s end after withdrawal

1

α m g β m i

2

α m g 2 m g β m i β m i 2

3

α m g 3 m g 2 β m i m g β m i 2 β m i 3

4

α m g 4 m g 3 β m i m g 2 β m i 2 m g β m i 3 β m i 4

5

α m g 5 m g 4 β m i m g 3 β m i 2 m g 2 β m i 3 m g β m i 4 β m i 5

We factor out β from every term except the first term.

For example, for Month 3, we have

f( 3 )=α m g 3 m g 2 β m i m g β m i 2 β m i 3 =α m g 3 β( m g 2 m i + m g m i 2 + m i 3 ) (5)

We highlight the exponents of mg and mi in every term. Table 3 displays the amount in the fund at the month’s end with highlighted exponents.

Table 3. The amount in the fund at the month’s end with highlighted exponents.

Month

The amount in the fund at the month’s end after withdrawal

1

α m g 1 β( m g 0 m i 1 )

2

α m g 2 β( m g 1 m i 1 + m g 0 m i 2 )

3

α m g 3 β( m g 2 m i 1 + m g 1 m i 2 + m g 0 m i 3 )

4

α m g 4 β( m g 3 m i + m g 2 m i 2 + m g 1 m i 3 + m g 0 m i 4 )

5

α m g 5 β( m g 4 m i + m g 3 m i 2 + m g 2 m i 3 + m g 1 m i 4 + m g 0 m i 5 )

2.4. Splitting Expressions for Every Month

Figure 2. Splitting of the expression (Example for Month 3).

For every month, the expression can be written as a subtraction of two distinct expressions: Expression 1 and Expression 2 (Figure 2).

2.5. Modeling Expression 1

For every month, Table 4 illustrates Expression 1 and the exponent of the variable mg in Expression 1.

Table 4. Expression 1 and the exponent of mg in Expression 1.

Month

Expression 1

Exponent of variable mg

1

a m g 1

1

2

a m g 2

2

3

a m g 3

3

4

a m g 4

4

5

a m g 5

5

Let the variable x represent the month number. For example, the value of x for Month 3 is 3.

Then, Expression 1 can be modeled by the function:

f( x )=α m g x (6)

2.6. Modeling Expression 2

For every month, it is evident that the expression that the variable β multiplies is a series (Figure 3).

Figure 3. The expression that β multiplies is a series (Example for Month 3).

Step 1: We isolate the series from the variable β.

Step 2: We isolate every term of the series.

Table 5. Isolation of every term of the series.

Month

Terms of the series

1st Term

2nd Term

3rd Term

4th Term

5th term

1

m g 0 m i 1

-

-

-

-

2

m g 1 m i 1

m g 0 m i 2

-

-

-

3

m g 2 m i 1

m g 1 m i 2

m g 0 m i 3

-

-

4

m g 3 m i 1

m g 2 m i 2

m g 1 m i 3

m g 0 m i 4

-

5

m g 4 m i 1

m g 3 m i 2

m g 2 m i 3

m g 1 m i 4

m g 0 m i 5

Step 3: We isolate the exponent of the variable mg and the exponent of the variable mi from every term of the series (Table 5).

Table 6. Isolation of the exponents of mg and mi from every term of the series.

Month

Term

Exponent of mg (j)

Exponent of mi (k)

Exponent of mg + Exponent of mi (j + k)

1

1st Term

0

1

1

2

1st Term

1

1

2

2nd Term

0

2

2

3

1st Term

2

1

3

2nd Term

1

2

3

3rd Term

0

3

3

4

1st Term

3

1

4

2nd term

2

2

4

3rd Term

1

3

4

4th Term

0

4

4

5

1st Term

4

1

5

2nd Term

3

2

5

3rd Term

2

3

5

4th Term

1

4

5

5th Term

0

5

5

Step 4: By observing Table 6, we deduce that the month number is equal to the number of terms in the series corresponding to that month. Then, the variable x (defined in Section 2.5) represents both the month number and the number of terms in the series corresponding to that month.

Step 5: Let j denote the exponent of the variable mg, and let k denote the exponent of the variable mi

Step 6: Therefore, by observing the trends in the values of j and k in Table 7, we formulate a generalized expression for representing the series for any arbitrary month x:

k=1 x m g j m i k

Step 7: By observing the rightmost column of Table 6, we can deduce that the sum of j and k for any arbitrary month x is equal to the month number x.

Then, we have:

j+k=x

j=xk (7)

We substitute j with x – k in the generalized expression for the series:

k=1 x m g xk m i k

Step 8: We combine the generalized expression for the series with the variable β, which was isolated in Step 1, to obtain the function modelling Expression 2:

f( x )=β( k=1 x m g xk m i k ) (8)

2.7. Formulating the Function

As shown in Figure 2, the main expression is formed by subtracting Expression 2 from Expression 1.

Therefore, the function modelling the main expression:

f( x )=a m g x β( k=1 x m g xk m i k ) (9)

We replace our simplifications mg and mi with their original expansions to obtain:

f( x )=α ( 1+ r g 100 ) x β( k=1 x ( 1+ r g 100 ) xk ( 1+ r i 100 ) k ) (10)

The function above is a function of the month number x, and f(x) represents the amount in the fund at the end of x months. The function models the fund of a systematic withdrawal plan, for which the variables α, β, rg, and ri can be inputted.

3. Stochastic Analysis and Simulation of Fund Longevity

3.1. Exploratory Data Analysis

The function that models the fund assumes that rg is constant. However, rg, which denotes the monthly growth rate of the S&P 500, is bound to vary. Therefore, we must account for the uncertainty in rg in the model. A dataset of 168 rg values was retrieved [4]. Each rg value corresponds to its respective month from January 2010 to January 2024. To check if the randomness in rg values follows a predictable pattern, exploratory data analysis is conducted by generating a histogram of rg values (Figure 4).

Figure 4. Histogram of rg values.

Observing the shape of the histogram, we can intuitively hypothesize that the variable rg is likely to be normally distributed. To validate our hypothesis, we shall check the normality of rg using a one-sample Kolmogorov–Smirnov test [5].

3.2. Kolmogorov–Smirnov Test

Hypotheses for the one-sample Kolmogorov–Smirnov test on rg:

H: rg is a random variable drawn from a normal distribution with μ equal to the mean of rg and σ equal to the standard deviation of rg.

HA: rg is not a random variable drawn from a normal distribution with μ equal to the mean of rg and σ equal to the standard deviation of rg.

Table 7. Kolmogorov-Smirnov test results.

Variable

Value

KS Statistic

0.0764

p-value

0.2582

μ

0.6933 (% per month)

σ

4.0266 (% per month)

The Kolmogorov–Smirnov test outputs a KS Statistic of 0.0764 and a p-value of 0.2582. Since the p-value of 0.2582 is more than the alpha value of 0.05, the null hypothesis is accepted. This comparison proves that rg is random variable normally distributed with μ as 0.6933 (% per month) and σ as 4.0266 (% per month).

3.3. Stochastic Simulations

Since we know that rg is a normally distributed random variable, we can run stochastic simulations for the function that models the fund. Each time we run a stochastic simulation; we note the number of months it takes for the function to fall below zero. In other words, we note how long the fund of a systematic withdrawal plan lasts. Let t be a new variable, where t refers to the number of months it takes for the fund to fall below zero. Since the function that models the fund depends on rg, the uncertainty in rg carries over to t. By running sufficient stochastic simulations, we can generate a histogram of t values to check if the randomness in t values follows a predictable pattern. We set realistic example values for α as $500000 and β as $5000 for the stochastic simulations. rg is sampled from a normal distribution with parameters (0.6933, 4.0266), and ri is 1.0021(Figure 5).

Figure 5. 10000 simulations: α = 500000, β = 5000, rg = normal (0.6933, 4.0266), ri = 1.0021.

Pseudocode for plotting a histogram of t values by running stochastic simulations:

loop i from 1 to 10000

f(x) = α

loop x while f(x) > 0

rg = np.random.normal(mu, sigma)

f(x) = (f(x)*(1 + (rg/100))) – (β*(1 + (ri/100))^x)

x++

end loop

t_values[i] = x

end loop

plot histogram of t_values

3.4. Outliers

Figure 6 suggests that outliers exist when t equals 10,000 months. These outliers are halted at 10,000 months due to the technical limitations of the R language. However, these outliers signify extreme cases when t approaches infinity. These cases happen when stochastic variations result in unrealistic consistently high rg values, which allow the fund to last infinitely. However, these cases are significantly rare, with a 0.58% (58 in 10,000 simulations) probability. Therefore, we remove the outliers to improve robustness, and we plot a new histogram.

Figure 6. Histogram of t values with outliers.

Figure 7. Histogram of t values without outliers.

The shape of the histogram of t values in Figure 7 suggests that the variable t is likely to be log-normally distributed. To validate our hypothesis, we shall log-transform [6] the t values and check the normality of the log-transformed t values using a one-sample Kolmogorov–Smirnov test.

Table 8. Kolmogorov-Smirnov test results.

Variable

Value

KS Statistic

0.0488

p-value

0.1318

μt

4.8466

σt

0.2721

The Kolmogorov–Smirnov test outputs a KS Statistic of 0.0488 and a p-value of 0.1318. Since the p-value is more than the alpha value of 0.05, the null hypothesis is accepted (Table 8). This comparison proves that t is a random variable drawn from a log-normal distribution with μt equal to 4.8466 and σt equal to 0.2721. Therefore, a continuous probability distribution function for t values, based on the log-normal distribution [7] is given by:

P( x )= 1 2 ( 1+erf( ln( x )4.8466 0.2721 2 ) ) (11)

Equation (11) applies specifically to the example given in this study. However, the log-normal nature of the distribution of t values holds for all values of α and β, as the shape of the distribution depends on rg, and rg follows a fixed normal distribution.

3.5. Percentiles

Table 9. Percentiles of t values.

Percentile

t Value (Months)a

10th

90

20th

101

30th

110

40th

119

50th

127

60th

136

70th

147

80th

160

90th

180

a. t values are rounded to the nearest integer.

Table 9 displays the percentiles of t values according to the log-normal distribution given by Equation (11). The deterministic model from Section 2, which does not account for randomness, predicts a t value of 137 months for α = 500,000, β = 5000, rg = 0.6933, and ri = 1.0021. rg was set to 0.6933, which is the mean and best estimate of rg, as the model requires a constant value. This prediction is validated by our stochastic analysis, as 137 is close to 127, which is the median of the log-normal distribution of t values. However, the deterministic model and the median are not a reasonable estimate of the fund’s longevity; they approximately provide a 50% safety margin, as 50% of the time, we will get a t value of less than 127 months, solely due to random chance. We must have a 90% safety margin, chosen in several equity funds, to have a high level of confidence in the fund’s longevity.

3.6. Stochastic Algorithm

Our goal is to develop an algorithm to output α, the initial fund requirement, with a 90% safety margin, based on inputs of β and tTarget Let tTarget be a new variable, which represents the number of months we require our fund to last.

Step 1: Given inputs β and tTarget, we use the deterministic model from Section 2 to find the 50% safety margin value for α.

When we set Equation (10) from Section 2 to equal zero, the x in the equation represents the number of months the fund of the systematic withdrawal plan lasts.

α ( 1+ r g 100 ) x β( k=1 x ( 1+ r g 100 ) xk ( 1+ r i 100 ) k )=0 (12)

We rearrange Equation (13) to obtain:

α ( 1+ r g 100 ) x =β( k=1 x ( 1+ r g 100 ) xk ( 1+ r i 100 ) k )

α= β( k=1 x ( 1+ r g 100 ) xk ( 1+ r i 100 ) k ) ( 1+ r g 100 ) x

α= β( k=1 t target ( 1+ r g 100 ) t target k ( 1+ r i 100 ) k ) ( 1+ r g 100 ) t target (13)

We plug the values of β, rg = 0.6933, ri = 1.0021, and tTarget into Equation (13) to calculate the 50% safety value of α.

Step 2: We run stochastic simulations to find the log-normal distribution of t values for the 50% safety value of α using the pseudocode from Section 3.3.

Step 3: We identify the 90th percentile t value from the log-normal distribution of t values, and we plug it into Equation (13) to find the 90% safety value of α.

Figure 8. Flowchart for the stochastic algorithm.

For example, we set the same example value from Section 3.3 of β as $5000, and we set the target t value as 127 months from Section 3.5. Running the stochastic algorithm from Figure 8 on these parameters returns α as $629,065. Compared to the lower initial α value of $500,000 from Section 3.3, which had a 50% safety margin, the higher α value of $629,065 suggested by our algorithm gives a 90% safety margin.

4. Conclusion

This study identified that the monthly growth rate of the S&P 500 is normally distributed, with a mean of 0.6933% per month and a standard deviation of 4.0266% per month. The algorithm performs efficiently for other stocks indexes like BSE SENSEX and FTSE 100, provided their monthly growth rate is normally distributed. Stochastic simulations using this normal distribution led to the observation of a log-normal distribution of the fund’s longevity after the removal of outliers. Based on the distribution of the fund’s longevity, this study proposes a stochastic algorithm that output the initial amount required to sustain a systematic withdrawal plan for a target duration with a 90% safety margin, given an initial monthly withdrawal amount. Applying the algorithm to countries with unstable inflation and stock fluctuations [8] may need larger safety margins of 95% to account for increased randomness. The algorithm outperforms conventional systematic withdrawal plan calculators since it is personalized by inputting the mean and standard deviation for any stock. A benefit of the algorithm is the removal of unrealistic favorable outliers. Therefore, the predictions of the initial fund requirements are not skewed by extreme data points [9]. Another advantage of the algorithm is its 90% safety margin, which ensures that retirement planners have a high level of confidence that the fund of their systematic withdrawal plan will last for their target durations.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1] Shah, T. and Baser, N. (2022) Global Mutual Fund Market: The Turn of the Month Effect and Investment Strategy. Journal of Asset Management, 23, 466-476.
https://doi.org/10.1057/s41260-022-00282-0
[2] Sodini, P. and Viceira, L.M. (2020) The Value of Diversification: Diversification is the Only Free Lunch in Finance. Harvard University.
https://scholar.harvard.edu/files/lviceira/files/ap7_annual_report-ps_and_lv-2020-01-29.pdf
[3] Macrotrends (2024) U.S. Inflation Rate 1960-2024. Macrotrends.
https://www.macrotrends.net/global-metrics/countries/USA/united-states/inflation-rate-cpi
[4] Investing.com (2024) S&P 500 Index Historical Data.
https://www.investing.com/indices/us-spx-500-historical-data
[5] Bristol University (2024) SPSS—Exploring Normality (Practical).
https://www.bristol.ac.uk/cmm/media/research/ba-teaching-ebooks/pdf/Normality%20-%20Practical.pdf
[6] Feng, C., Wang, H., Lu, N., Chen, T., He, H., Lu, Y. and Tu, X.M. (2014) Log-Transformation and Its Implications for Data Analysis. Shanghai Archives of Psychiatry, 26, 105-109.
https://doi.org/10.3969/j.issn.1002-0829.2014.02.009
[7] StatProofBook (2024) Proof: Cumulative Distribution Function of the Log-Normal Distribution.
https://statproofbook.github.io/P/lognorm-cdf.html
[8] The Global Economy (2021) Stock Price Volatility—Country Rankings.
https://www.theglobaleconomy.com/rankings/Stock_price_volatility/
[9] Sullivan, J.H., Warkentin, M. and Wallace, L. (2021) So Many Ways for Assessing Outliers: What Really Works and Does It Matter? Journal of Business Research, 132, 530-543.
https://doi.org/10.1016/j.jbusres.2021.03.066

Copyright © 2025 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.