This post will introduce the first part (of multiple) where we build up a personal finance model to help simulate future time periods based on certain chosen input variables. We will input variables such as our current investable asset base, our annual salary, expected monthly inflows and outflows and a range of other relevant values. Firstly, after our necessary imports, we look to start on modelling inflows over time.
import pandas as pd import numpy as np import random import pandas_datareader.data as web import datetime import matplotlib as mpl import matplotlib.pyplot as plt %matplotlib inline plt.rcParams['figure.figsize'] = (12,6) # Set up model inflows inflows = {'active_annual_income':50_000} variables = {'start_date' : "01/01/2020", 'years': 10} income_gains_storage = [] months = variables['years'] * 12 for month in range(months): income = inflows['active_annual_income'] / 12 income_gains_storage.append(income) plt.plot(pd.Series(income_gains_storage).cumsum()) plt.show()

This looks like a pretty simplistic results, the cumulative income we have earned over 10 years, at 50k annual salary is unsurprisingly 10 x 50,000 = 500,000. One thing we might want to introduce is some logic to account for the fact the nasty tax man gets his hands on at least some (way too much) of it. For the moment, for simplicity’s sake, I will assume a flat tax rate of 25%. Of course most countries have a sliding scale of income bands which are subject to increasing rates of tax the more one earns. We can come back to this later and see if we can implement something more along those lines – but for now lets stick with our flat 25%.
Let’s add an entry to our ‘variables’ dictionary to represent the tax rate.
I have purposefully named it to reflect the fact this is tax being applied to one’s income from employment – what i have called “active_annual_income” in the “inflows” dictionary. This is as opposed to taxes that we might incur against other sources of income, such as investment income (which we shall introduce further down the line).
inflows = {'active_annual_income':50_000} variables = {'start_date' : "01/01/2020", 'years': 10, 'tax_on_active_income_gains': 0.25} # add our tax rate on active income in % terms income_gains_storage = [] months = variables['years'] * 12 for month in range(months): # add the effect of the tax rate being applied to our income income = (inflows['active_annual_income'] * (1 - variables['tax_on_active_income_gains'])) / 12 income_gains_storage.append(income) plt.plot(pd.Series(income_gains_storage).cumsum()) plt.xlabel('Month') plt.ylabel('Cumulative Income') plt.show()

So again, at this point, nothing too fantastically insightful or useful in terms of output; we just end up now after the 10 years with 375K, which is just the 75% left over after the taxman has taken his 25%. This logic is, of course, very over-simplistic as not only do we not account for any kind of sliding-scale tax bands at which different levels of income are taxed, but we also have no concept of the time-value-of-money here. A dollar of income taken by the tax man today hurts us more than a dollar taken in say 5 years from now. Not only would we be able to earn interest on the dollar if we were allowed to keep it, we could also invest it and put it to work in search of higher returns.
We will get around to modelling these investment returns shortly – in the real-world they are stochastic in nature and display a certain element of random variation from one period to the next. We can’t just assume if our calculated annual expected investment returns are 10% per year that we will realise those exact returns every year, nor receive them in a linear fashion of 12 nice, evenly sized monthly instalments. The fact is that it is not only the size of our returns that determine our ending value of wealth accumulated; the order of those returns matter, along with their volatility. I won’t spend too much time going through geometric means vs arithmetic means and such, but suffice to say the higher the volatility of a set of returns, the more damaging that is to our end wealth (on average). As an example if we were to see returns of -20% in period 1 and +20% in period 2 we would end up with 0.8 x 1.2 = 0.96 i.e. we would be down 4% from a starting capital of 1, vs if we saw more volatile returns of -50% and +50 where we would end up with 0.5 x 1.5 = 0.75 – i.e. we would be down 25%. This is even though both examples give us an arithmetic return of 0% (-0.2 + 0.2 = 0) and (-0.5 + 0.5 = 0). If interested in further reading, a more detailed explanation of the relationship between arithmetic and geometric means can be found here: https://brilliant.org/wiki/arithmetic-mean-geometric-mean/
Just before we move on to the investment logic, let’s quickly add some logic to allow for an annual salary increase, set as a percentage value of the current active annual salary. The way we shall implement it is to assume our boss will sit down with us once a year and offer us an increase that will take effect at that time. This is how we go about it:
inflows = {'active_annual_income':50_000} variables = {'start_date' : "01/01/2020", 'years': 10, 'tax_on_active_income_gains': 0.25, 'avg_ann_income_raise':0.05} # add our annual salary raise in % terms income_gains_storage = [] months = variables['years'] * 12 for month in range(months): income = (inflows['active_annual_income'] * (1 - variables['tax_on_active_income_gains'])) / 12 income_gains_storage.append(income) # every 12 months we increase our base salary by the average annual percentage increase if (month % 12 == 0) and (month > 0): # we dont want to apply the raise on the very first month inflows['active_annual_income'] *= (1 + (variables['avg_ann_income_raise'])) plt.plot(pd.Series(income_gains_storage).cumsum()) plt.xlabel('Month') plt.ylabel('Cumulative Income') plt.show()

Ok now let’s get to modelling those stochastic investment returns mentioned previously. There are numerous ways to go about this task; the exact approach may well differ based on individual circumstances or beliefs etc. To keep it relatively simple for illustration purposes, we shall use the S&P 500 Total Return index as a proxy for “the market” and use historic data to calculate our values regarding “average annual returns” and related levels of volatility. We will calculate these two values (average return and standard deviation) – we shall then use them as inputs to draw from a Standard Normal distribution and use the values received as our assumed “market returns” across the period in question.
So we are now going to:
- Download historic S&P 500 index prices
- Calculate the average monthly return and standard deviation of a chosen historic period
- Store these values to use later as inputs to our Standard Normal distribution
The values are coming out as a monthly return of 0.6% and a monthly volatility (standard deviation) of 4.15%.
Before we attack this logic, we need to create a new input in our “inflows” dictionary to represent our starting capital, or starting assets. Also we add the calculated values of average monthly return and volatility to our “variables” dictionary.
# download S&P historic prices data start, end = datetime.datetime(2000, 12, 31), datetime.datetime(2020, 1,1) tickers = ["^SP500TR"] sp = pd.DataFrame([web.DataReader(ticker, 'yahoo', start, end).loc[:, 'Adj Close'] for ticker in tickers], index=tickers).T.fillna(method='ffill') # calculate average monthly return & volatility sp_monthly_pct_return = sp.resample('M').last().pct_change().mean().values[0] sp_monthly_std_dev = sp.resample('M').last().pct_change().std().values[0] inflows = {'active_annual_income':50_000, 'starting_assets': 250_000} # add value of starting assets variables = {'start_date' : "01/01/2020", 'years': 10, 'tax_on_active_income_gains': 0.25, 'avg_ann_income_raise':0.05, 'avg_monthly_market_returns': sp_monthly_pct_return, # add market return data to our inputs 'avg_monthly_market_volatility': sp_monthly_std_dev} # add market volatility data to our inputs income_gains_storage = [] investment_gains_storage = [] # create list to store our investment returns # create a list to store each period's starting assets value assets_starting_list = [inflows['starting_assets']] # input the first value in the list as today's assets assets_ending_list = [] # create a list to store each period's ending assets value months = variables['years'] * 12 for month in range(months): # small bit of logic to see if this is our first time through the monthly loop # and if it is not, then start to use the ending assets from the previous # period as this period's starting assets. if assets_ending_list: assets_starting_list.append(assets_ending_list[-1]) income = (inflows['active_annual_income'] * (1 - variables['tax_on_active_income_gains'])) / 12 income_gains_storage.append(income) if (month % 12 == 0) and (month > 0): inflows['active_annual_income'] *= (1 + (variables['avg_ann_income_raise'])) # generate a monthly market return by drawing from a normal distribution market_return = np.random.normal(variables['avg_monthly_market_returns'], variables['avg_monthly_market_volatility'], 1)[0] # calculate the investment return investment_return = assets_starting_list[-1] * market_return # store the investment return value in a list investment_gains_storage.append(investment_return) # calculate the end of period assets value assets_ending = assets_starting_list[-1] + investment_return + income # store ending assets value assets_ending_list.append(assets_ending) plt.plot(pd.Series(investment_gains_storage).cumsum()) plt.xlabel('Month') plt.ylabel('Cumulative Investment Returns') plt.show()

Now when we plot our series of investment returns we can see that they are indeed stochastic in nature, with the inherent randomness of returns meaning some months we are up, some we are down even though our average monthly return used as an input to the Standard Normal function was positive.
As a slight point to note – please be aware that some of the following plots are showing the series of cumulative monthly investment returns, whereas some of them are showing the series of ending period asset value – they are different things! Apologies if that is obvious, just I appreciate the plots can often look very similar on first sight. Make sure you check the y-axis label on each plot.
If you run the above code a number of times, each outcome will be different due to that fact – this will lead us on nicely later to the world of Monte Carlo methods and creating distributions of ending values – I will leave the juicy details until we get around to that section.
Ok so now we have modelled our monthly active salary (with an annual percentage increase), along with our monthly investment returns which we have subjected to a stochastic element to mimic the real-world.
If we plot our series of ending asset values and note how much we have ended up with, you may assume we are being over-optimistic. You would, of course, be correct; firstly we haven’t applied any taxes to our investment gains, and also the elephant in the room is that we are yet to include any kind of “outflows” – they’re those pesky bills and living costs we have to shell out to do things like eat food and have somewhere to live. Once we start to apply the “debit” side of the equation, we will see the situation balance out somewhat (and here’s me thinking I’d be a multi-millionaire in a few years).
plt.plot(pd.Series(assets_ending_list)) plt.xlabel('Month') plt.ylabel('Ending Asset Value') plt.show()

So let’s first apply a tax rate to our investment returns – the process is similar to when we applied tax to our active income. First we add an investment tax rate to the variables dictionary – the code is shown below:
start, end = datetime.datetime(2000, 12, 31), datetime.datetime(2020, 1,1) tickers = ["^SP500TR"] sp = pd.DataFrame([web.DataReader(ticker, 'yahoo', start, end).loc[:, 'Adj Close'] for ticker in tickers], index=tickers).T.fillna(method='ffill') sp_monthly_pct_return = sp.resample('M').last().pct_change().mean().values[0] sp_monthly_std_dev = sp.resample('M').last().pct_change().std().values[0] inflows = {'active_annual_income':50_000, 'starting_assets': 250_000} variables = {'start_date' : "01/01/2020", 'years': 10, 'tax_on_active_income_gains': 0.25, 'avg_ann_income_raise':0.05, 'tax_on_investment_gains': 0.35, # add investment returns tax rate 'avg_monthly_market_returns': sp_monthly_pct_return, 'avg_monthly_market_volatility': sp_monthly_std_dev} income_gains_storage = [] investment_gains_storage = [] assets_starting_list = [inflows['starting_assets']] assets_ending_list = [] months = variables['years'] * 12 for month in range(months): if assets_ending_list: assets_starting_list.append(assets_ending_list[-1]) income = (inflows['active_annual_income'] * (1 - variables['tax_on_active_income_gains'])) / 12 income_gains_storage.append(income) if (month % 12 == 0) and (month > 0): inflows['active_annual_income'] *= (1 + (variables['avg_ann_income_raise'])) market_return = np.random.normal(variables['avg_monthly_market_returns'], variables['avg_monthly_market_volatility'], 1)[0] # apply the tax rate to investment gains investment_return = (assets_starting_list[-1] * market_return) * (1 - variables['tax_on_investment_gains']) # store the investment return value in a list investment_gains_storage.append(investment_return) # calculate the end of period assets value assets_ending = assets_starting_list[-1] + investment_return + income # store ending assets value assets_ending_list.append(assets_ending) plt.plot(pd.Series(investment_gains_storage).cumsum()) plt.xlabel('Month') plt.ylabel('Cumulative Investment Returns') plt.show()

Let’s now start to add some regular monthly outflows and living costs to even up the balance. I have gone about it by creating a new “outflows” dictionary as shown below. I have also introduced a variable to reflect the general level of inflation which we will apply to our costs on an annual basis (much like we hiked our salary by a certain percentage once a year). We could of course deal with the annual inflation level in a stochastic manner as we did the investment returns, but for simplicity here we shall just apply them deterministically.
start, end = datetime.datetime(2000, 12, 31), datetime.datetime(2020, 1,1) tickers = ["^SP500TR"] sp = pd.DataFrame([web.DataReader(ticker, 'yahoo', start, end).loc[:, 'Adj Close'] for ticker in tickers], index=tickers).T.fillna(method='ffill') sp_monthly_pct_return = sp.resample('M').last().pct_change().mean().values[0] sp_monthly_std_dev = sp.resample('M').last().pct_change().std().values[0] inflows = {'active_annual_income':50_000, 'starting_assets': 250_000} # add our outflows dictionary outflows = {'rent':1500, 'credit_card_payment':750, 'medical_insurance':250, 'pension_contribution':500, 'misc': 1500} variables = {'start_date' : "01/01/2020", 'years': 10, 'tax_on_active_income_gains': 0.25, 'avg_ann_income_raise':0.05, 'avg_ann_inflation': 0.02, # add annual inflation rate 'tax_on_investment_gains': 0.35, 'avg_monthly_market_returns': sp_monthly_pct_return, 'avg_monthly_market_volatility': sp_monthly_std_dev} income_gains_storage = [] investment_gains_storage = [] assets_starting_list = [inflows['starting_assets']] assets_ending_list = [] months = variables['years'] * 12 for month in range(months): if assets_ending_list: assets_starting_list.append(assets_ending_list[-1]) # create a variable to keep track of our assets value assets = assets_starting_list[-1] # first we subtract our monthly costs from our overall assets assets -= (outflows['rent'] + outflows['credit_card_payment'] + \ outflows['medical_insurance'] + outflows['pension_contribution'] + \ outflows['misc']) # move our calculation of investment returns above that for income # and use our new "assets" variable as the investment base market_return = np.random.normal(variables['avg_monthly_market_returns'], variables['avg_monthly_market_volatility'], 1)[0] investment_return = (assets * market_return) * (1 - variables['tax_on_investment_gains']) investment_gains_storage.append(investment_return) # add investment returns to our "assets" variable assets += investment_return income = (inflows['active_annual_income'] * (1 - variables['tax_on_active_income_gains'])) / 12 income_gains_storage.append(income) if (month % 12 == 0): inflows['active_annual_income'] *= (1 + (variables['avg_ann_income_raise'])) # increase our outflows by inflation each year # I have chosen to list each one individually to enable easy customisation # (rather than looping over the outflows dict to increase them) outflows['rent'] *= (1 + (variables['avg_ann_inflation'])) outflows['credit_card_payment'] *= (1 + (variables['avg_ann_inflation'])) outflows['medical_insurance'] *= (1 + (variables['avg_ann_inflation'])) outflows['pension_contribution'] *= (1 + (variables['avg_ann_inflation'])) outflows['misc'] *= (1 + (variables['avg_ann_inflation'])) # add income gain to our "assets" variable assets += income # calculate the end of period assets value assets_ending = assets # store ending assets value assets_ending_list.append(assets_ending) plt.plot(pd.Series(assets_ending_list)) plt.xlabel('Month') plt.ylabel('Ending Asset Value') plt.show()

Above I have moved the application of outflows to the asset base above the code that adds the income or calculates the investment returns. The thinking here is that if you have bills to pay, then you make sure you pay them first as soon as you can, or at least hold the money aside for the month so you are certain you will have it when needed – rather than invest it and risk not having enough when it comes time to pay. I also assume here that salaries will be paid at the end of the month – or at least not early enough to be sat in your bank account, waiting to be invested or able to be used to pay bills for that month (better safe than sorry).
There is one more bit of logic I would like to apply before we move forward; currently the code allows for a situation where we go “bust”, losing everything we have in terms of investable assets, and continue “investing” a negative amount. This is obviously nonsensical as you cant make positive returns from a negative asset base, just because the market returns are negative. We will add a flag that checks to make sure our asset base is in positive territory before we start calculating and applying investment returns; if our asset base is negative, there are a few things we can do to ac count for that – just for example we could:
- Consider that situation “unnaceptable” to the extent we class ourselves as “ruined/bankrupt” and end the simulation with a “failure” flaf or label.
- Allow the simulation to continue and collect salary income each month, but each month where assets are non-positive, we set the investment gains to zero by default.
- Apply some non-trivial logic that applies interest and penalty charges, allow us to borrow and re-invest etc.
The last option, in my view, is just not realistic nor really the point of this exercise, with the best options being between 1 and 2. Which ever you choose is up to you – I tend to go for option 1, considering an asset base of zero to be unacceptable and a sign of abject failure with regards to money management and such. You could even argue that even getting close to a zero value asset base is unacceptable and want to set the “failure” threshold somewhere above zero. We can revisit these options later perhaps if it is of benefit.
I deal with a failed simulation by recording a “failure flag” – mark that run as failed for later analysis. This reasoning will become a bit clearer later on when we start to apply Monte Carlo methods to our model, carrying out 100s, 1000s or more simulations per “experiment”. Don’t worry if what I am saying isn’t making sense – we will tackle the Monte Carlo element and explain in due course.
Let’s apply the logic now to stop the simulation if we lose our entire asset base and effectively become financially “ruined”.
start, end = datetime.datetime(2000, 12, 31), datetime.datetime(2020, 1,1) tickers = ["^SP500TR"] sp = pd.DataFrame([web.DataReader(ticker, 'yahoo', start, end).loc[:, 'Adj Close'] for ticker in tickers], index=tickers).T.fillna(method='ffill') sp_monthly_pct_return = sp.resample('M').last().pct_change().mean().values[0] sp_monthly_std_dev = sp.resample('M').last().pct_change().std().values[0] inflows = {'active_annual_income':50_000, 'starting_assets': 75_000} # add our outflows dictionary outflows = {'rent':1500, 'credit_card_payment':750, 'medical_insurance':1250, 'pension_contribution':500, 'misc': 1500} variables = {'start_date' : "01/01/2020", 'years': 10, 'tax_on_active_income_gains': 0.25, 'avg_ann_income_raise':0.05, 'avg_ann_inflation': 0.02, 'tax_on_investment_gains': 0.35, 'avg_monthly_market_returns': sp_monthly_pct_return, 'avg_monthly_market_volatility': sp_monthly_std_dev} income_gains_storage = [] investment_gains_storage = [] assets_starting_list = [inflows['starting_assets']] assets_ending_list = [] months = variables['years'] * 12 # create a variable to signal if we have become financially "ruined" ruined = False for month in range(months): if assets_ending_list: assets_starting_list.append(assets_ending_list[-1]) assets = assets_starting_list[-1] assets -= (outflows['rent'] + outflows['credit_card_payment'] + \ outflows['medical_insurance'] + outflows['pension_contribution'] + \ outflows['misc']) # check if our asset base has a positive value and, if # not then set the "ruined" flag to 1 and end the simulation if assets <= 0: inv_gain = 0 ruined = True break market_return = np.random.normal(variables['avg_monthly_market_returns'], variables['avg_monthly_market_volatility'], 1)[0] investment_return = (assets * market_return) * (1 - variables['tax_on_investment_gains']) investment_gains_storage.append(investment_return) assets += investment_return income = (inflows['active_annual_income'] * (1 - variables['tax_on_active_income_gains'])) / 12 income_gains_storage.append(income) if (month % 12 == 0): inflows['active_annual_income'] *= (1 + (variables['avg_ann_income_raise'])) outflows['rent'] *= (1 + (variables['avg_ann_inflation'])) outflows['credit_card_payment'] *= (1 + (variables['avg_ann_inflation'])) outflows['medical_insurance'] *= (1 + (variables['avg_ann_inflation'])) outflows['pension_contribution'] *= (1 + (variables['avg_ann_inflation'])) outflows['misc'] *= (1 + (variables['avg_ann_inflation'])) assets += income assets_ending = assets assets_ending_list.append(assets_ending) plt.plot(pd.Series(assets_ending_list)) plt.xlabel('Month') plt.ylabel('Ending Asset Value') plt.show()

I have changed the values in the inflow and outflow dictionaries above to show a simulation that ends in “ruin” before the whole 10 year period has been simulated (I reduced the starting_assets in the inflow dictionary, and increased the medical insurance payment in the outflows dictionary FYI). We can see the x-axis only reaches as far as 36 or so months rather than the full 120 month (10 year) simulation period stated in the variables dictionary.
We will record and store the “ruined” flag a bit later on, in fact we will revisit this in the next post as this is getting a bit long as it is. Part 2 will cover a variety of things, e.g. Monte Carlo methods and subsequent probability estimates, intriducing “scenarios” with various “shock factors” implements, e.g. a one off asset-base reduction/loss, unexpected job and income loss, increase in predicted inflation costs among others.
Until next time…
9 comments
If ‘active_annual_income’:50_000 and ‘tax_on_active_income_gains’: 0.25 then net_income per month is 3_125.
If outflows = {‘rent’:1500,
‘credit_card_payment’:750,
‘medical_insurance’:1250,
‘pension_contribution’:500,
‘misc’: 1500}
Then total outflows per month is 4_500.
That means net income per month do not cover outflos.
Graph cannot be like showed.
You are missing the investment gains on the starting asset base of 250,000 in the charts which end in profit. You can see in the last chart when you only have a starting asset base of 75,000 you do indeed end in “ruin”, i.e. loss as the inflows do not cover the outflows. The charts are correct.
If outflows is 1500 per month bigger than income then assets needs to sell to cover that difference.
Or one can try to get a loan for that. That increases monthly outflows, and so on…
You are misunderstanding – you receive investment income on your asset base. Your monthly income is not 3_125 as you think; you need to add to that the investment income you make each month.
From that investment income one must use 1375 for living costs because costs are 4500 and salary covers only 3125 of the living costs. 1375 from 250 000 assets means that gain must be at least about 6.6% which is quite big percentage for normal investment. So propably assets gets smaller every month because gain is not enough big.
Again I have to disagree, I have modelled the investment returns stochastically using the actual volatility and return characteristic s of the S&P over the period starting in the year 2000. 6% is well within normal bounds and historic norms for annual equity returns. I take your point that nothing is guaranteed, hence the very reason we are running Monte Carlo simulations to get an idea of the range of possible outcomes and eventualities.
I have also modelled it such that outflows are subtracted before any investment returns are calculated as you are correct, you can’t earn returns on money that you have paid out on bills etc.
Mostly I get results like this:
starting_assets investment_gains income outflow ending_assets
0 $250,000.00 $-6,682.30 $3,125.00 $4,500.00 $241,942.70
1 $241,942.70 $-10,934.31 $3,125.00 $4,500.00 $229,633.39
2 $229,633.39 $-1,059.64 $3,125.00 $4,500.00 $227,198.74
3 $227,198.74 $3,181.61 $3,125.00 $4,500.00 $229,005.35
4 $229,005.35 $-1,409.42 $3,125.00 $4,500.00 $226,220.93
…
115 $120,169.61 $-2,645.99 $4,847.90 $5,377.92 $116,993.60
116 $116,993.60 $529.42 $4,847.90 $5,377.92 $116,993.00
117 $116,993.00 $6,162.87 $4,847.90 $5,377.92 $122,625.86
118 $122,625.86 $-2,676.01 $4,847.90 $5,377.92 $119,419.84
119 $119,419.84 $-33.18 $4,847.90 $5,377.92 $118,856.64
There is also some better results, but about eight runs from ten gives negative results.
So two of the ten are winning in the ‘game’.
Well, that’s normal in investment ‘game’.
Well now we are having a different conversation – I am not suggesting that one SHOULD build up monthly bills and outflows of that amount just because they have an asset base of 250k. But the point remains that the chart is indeed generated from the scenario where those values are the inputs. I have just rerun the code a few times and it only took a couple of runs to confirm that.
As a test I run the program 100 times and I got following results
Assets_ending
count 100
mean $236,361.11
std $98,095.49
min $41,482.41
25% $167,570.80
50% $226,733.57
75% $282,030.96
max $577,063.04
Testing several times I found that mean is about 250 000 and media below 250 000.
Minimum can be low and maximum can be high.