Solve for withdrawal rate in Monte Carlo simulation of retirement I've been working with compound returns and distribution of wealth over time for quite some time now and I feel like I am hitting a wall.
What am I trying to achieve? Imagine that you are about to retire and have saved $100 000. You are planning to withdraw your money monthly over the next 20 years. If we assume that returns are normally distributed with:
$${\mu} = 7 \%$$
$${\sigma} = 15 \%$$
I can simulate $1000$ returns from the above distribution in month $1$. Subtract $X$ from my initial pot, realise a return and calculate the new value of the pot at the end of that month. This gives me $1000$ values at the end of period 1 and I can take the median of all these. Next month the process is repeated. I simulate 1000 returns for month 2 now, subtract the same amount $X$ from the opening balance (end of period 1 values) and multiply by the returns. This gives me yet another 1000 values at the end of period 2 and I can take the median of those. I keep doing this for all the subsequent months and this gives me the path of all median values over $20$ years with withdrawals. There has to exist a way to calculate the amount X (which is constant in each period) that allows me to exhaust my retirement pot (terminal value of 0) at the end of the withdrawal period.
From Farago, Adam and Hjalmarsson, Erik, Compound Returns we know that compound returns have the following properties:
$${\psi} = \log\left(\frac{\mu^2}{\sqrt{\sigma^2+\mu^2}}\right)%$$
and
$${\eta} = \sqrt{\log\left(\frac{\sigma^2}{\mu^2}+1\right)} %$$
and we can find the $\alpha$-quantile of compound returns with:
$$q_\alpha\left(X_T\right) = e^{T\psi+\sqrt{T}\eta\Phi^{-1}\left(\alpha\right)}$$
i.e. I can use this formula to get the 25th percentile of compound returns from the above distribution in period T = 10
Unfortunately, this formula does not hold when you introduce either monthly savings or monthly withdrawals to the picture, as the distribution of compound wealth is then changed.
I'd appreciate if anyone can give me any pointers on solving for withdrawal amount X as per the above.
 A: It seems your problem can be easily solved by systematic trial and error. But your problem is not fully specified yet. Your terminal wealth is a random variable, i.e. it will depend on the performance of your investment in each scenario. This is why you need to define an additional criterion. Examples could be zero expected ultimate value, or caps on the shortfall probability or the application of a utility function.
Let us assume you have defined this criterion, which means for each terminal wealth you know whether the result is "too good" i.e. the withdrawal amount was too small, or too bad, i.e. the amount was too large. Then the trial and error solution works as follows:

*

*Define a reasonable initial guess for the withdrawal amount $x_0.$

*Simulate the 20 years with withdrawal $X=x_0$ a couple of times (a few hundred or few thousand times depending on your criterion).

*Look at the empirical distribution of your criterion. Decide whether the average of your criterion is too good, then increase $X$ or too bad then decrease the withdrawal.

*Repeat 2. and 3. until you are satisfied.

This should be easily doable for all "reasonable" criteria whose expectation is monotonic in the withdrawal rate.
Two further remarks:

*

*I do not see why calculating the median is useful.

*7% expected return???

Distribution of ultimate wealth but no closed form
It is possible to obtain an expression for the distribution of ultimate wealth. This shows that no closed form solution exists. Denote the initial wealth by $W$ the withdrawal amount by $x$. Let $t=1,\ldots,T$ be the periods until ultimate and $r_t$ the period (i.e. monthly) return. Now define the accumulated account value $A_t$ as
$$ A_{t+1}=(1 + r_{t+1}) A_t\quad\text{ with }A_0=1.$$
Your final wealth without withdrawals is
$$ W_T=W A_T$$
and the final value of x units withdrawal at time $t$ is
$$ X_t = x \frac{A_T}{A_t}.$$
So the distribution of final wealth is simply
$$ W_T-\sum_{t=1}^T X_t=A_T(W-x\sum_{t=1}^T\frac{1}{A_t}).$$
The variable $A_t$ is lognormal and it is well known that sums of lognormals do not have closed form expressions. Although various analytic expressions do exist see e.g. this post or that.
