Python Monte Carlo vs Bootstrapping

In this article I thought I would take a look at and compare the concepts of “Monte Carlo analysis” and “Bootstrapping” in relation to simulating returns series and generating corresponding confidence intervals as to a portfolio’s potential risks and rewards.

Both methods are used to generate simulated price paths for a given asset, or portfolio of assets but they use slightly differing methods, which can appear reasonably subtle to those who haven’t come across them before. Technically Bootstrapping is a special case of the Monte Carlo simulation, hence why it may seem a little confusing at first glance.

With Monte Carlo analysis (and here we are talking specifically about the “Parametric” Monte Carlo approach) the idea is to generate data based upon some underlying model characteristics. So for example, we generate data based upon a Normal distribution, specifying our desired inputs to the model, in this case being the mean and the standard deviation. Where do we get these input figures from I hear you ask…well more often than not people tend to use values based on the historic, realised values for the assets in question.

So if we were attempting to run some parametric Monte Carlo runs to generate simulated data for say Apply stock, we would tend to measure and calculate the mean and standard deviation of the stocks actual historic returns over a length of time, and use those as inputs to the model. This is one of the weaknesses of this approach, in as much as the model output and corresponding inferences are reliant on the assumption that future returns will display the same characteristics of historic returns (at least those that were used to calculated the model inputs).

So what is Bootstrapping and how does it differ? Well Bootstrapping also uses historic returns as a model input, but this time they are used more explicitly. Rather than calculating the underlying characteristics of the returns and then plugging those into a parametric model, we actually generate our data by sampling from the historic return distribution itself.

It is important to note here that Bootstrapping in involves “replacement” and falls under the notion of a “sampling with replacement” method.

That means that when a random sample is extracted from the historic return distribution, it is not “thrown away” and removed from the “hat” as it were, but rather it is replaced and put back in order that it may possibly be chosen again during the following sampling extractions.

This is a crucial point to note as it results in a fundamentally different outcome that if one was to “sample without replacement” and each data point were to be removed from the sample once chosen at any point.

The logic behind the Bootstrapping method is that if we use sampling with replacement, then each sample that is drawn, if random, will have the same chance of appearing as it would in “real life” – i.e. as it would in the actual markets for that particular stock (this again relies on the assumption that the future return distribution will retain the same characteristics as the historic return distribution the samples are being drawn from. That is to say both the future and past return distributions are drawn from the same “population”).

So that explanation out of the way, I thought I would collect data for a basket of of assets, create an equally weighted portfolio and then run both parametric Monte Carlo and Bootstrapping simulations, then compare the two results – see how similar of different they end up being!

First some basic imports and miscellaneous variable settings (stylesheet for charts, “magic” call to enable inline plotting for matplotlib in Jupyter notebooks, and setting the “figsize” variable for use in my charting calls).

import pandas as pd
import numpy as np
from functools import reduce
import pandas_datareader.data as web
import datetime
import random
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns
%matplotlib inline
mpl.style.use('ggplot')
figsize = (15, 8)

Now we download some price data from Yahoo for various stock indices using “pandas_datareader”, and rebase them all to start at 1 for purposes of comparison.

start, end = datetime.datetime(2009, 12, 30), datetime.datetime(2019, 5, 29)
tickers = ["^DJI", "^IXIC", "^GSPC", "^STOXX50E", "^N225", "^GDAXI"]
asset_universe = pd.DataFrame([web.DataReader(ticker, 'yahoo', start, 
                     end).loc[:, 'Adj Close'] for ticker in tickers],
                     index=tickers).T.fillna(method='ffill')
asset_universe = asset_universe/asset_universe.iloc[0, :]

Let’s plot our price series for our downloaded assets.

asset_universe.plot(figsize=figsize)

Now when carrying out Bootstrapping on a portfolio of assets, it is vitally important to make sure that we proceed correctly. Our approach must account for any correlations between the assets, as if we don’t we will get results that deviate from reality.

Take for example two stocks that are very strongly negatively correlated, if we were to sample independently for each stock when making a random draw, we may draw a sample that occurred on one particular day for stock 1, and a draw for stock 2 that occurred on a different day. Well if we do this, and we were to see that our draw for stock 1 was strongly positive, as was that for stock 2 could we really trust this to be a situation that was truly representative of the real relationship between the stocks?

The answer is no, because we are comparing “apples with oranges” as it were – we need the samples for our portfolio constituents to be drawn from the same time period, only then will the correlations inherent between all the assets be correctly captured in our random samples.

So we could proceed by generating multiple random draws (with replacement) from all our portfolio constituents individual historic return series then weighting them accordingly, and then finally summing the weighted returns and recording the corresponding output as our Bootstrapped “portfolio return”. We would then repeat this procedure many times, each time recording the simulated “portfolio return” – this collection of simulated return paths would be our Bootstrapped output.

Alternatively, we could build the portfolio return by weighting the constituent historic returns accordingly, summing them up and then carry out the Bootstrapping process on that single portfolio historic return distribution. The results are pretty much analogous as the way the portfolio was constructed in the second approach also inherently retains the effect of any correlations between the constituent assets as the return series was calculated using weighted constituent returns that happened on the same day. So we could just Bootstrap that single portfolio, again generating multiple simulated return paths and the collection of those paths would be our Bootstrapped output.

Let’s begin with the second approach and create our equally weighted portfolio return series. We just take the mean of the individual constituent returns for an equally weighted portfolio – simple as that. Let’s then plot the “price series” of our portfolio against the individual constituents.

portfolio_returns = asset_universe.pct_change().dropna().mean(axis=1)
portfolio = (asset_universe.pct_change().dropna().mean(axis=1) + 1).cumprod()
asset_universe.plot(figsize=figsize, alpha=0.4)
portfolio.plot(label='Portfolio', color='black')
plt.legend()
plt.show()

As you would expect the returns end up being somewhere in the middle of the individual returns. In fact as this is an equally weighted portfolio it will by definition end up exactly in the “middle” of the constituent returns.

Let’s now carry out the Bootstrapping process on our portfolio return series and plot the results.

portfolio_bootstrapping = (1+pd.DataFrame([random.choices(list(
    portfolio_returns.values), k=252) for i in 
    range(1000)]).T.shift(1).fillna(0)).cumprod()
portfolio_bootstrapping.plot(figsize=figsize, legend=False, linewidth=1, alpha=0.2, color='b')

And just to prove that approaching it from the other way around is analogous we do that below. We take samples of our individual constituent returns series and use them to create our Bootstrapped simulations. We end with the same result (there or there abouts, of course there is a random element which will make every simulation different even if based on the same approach) – below is the code that achieves that.

asset_universe_returns = asset_universe.pct_change()
portfolio_constituents_bootstrapping = 
    pd.DataFrame([((asset_universe_returns.iloc[random.choices(
    range(len(asset_universe)), k=252)]).mean(axis=1)+1).cumprod().values 
    for x in range(1000)]).T
portfolio_constituents_bootstrapping.plot(figsize=figsize, legend=False, linewidth=1, alpha=0.2, color='purple')

Finally we look to use the parametric Monte Carlo method, after which we can run a quick comparison of the results between the various approaches.

AS stated previously, the parametric Monte Carlo method involves using the characteristics of the underlying population to generate random samples of values. The characteristics we are talking about here are the mean and standard deviation (or variance) of the historic return distribution. These values will then be fed into a model that randomly samples from a normal distribution with mean and standard deviation equal to that of the historic returns.

Let’s extract those figures first for our combined portfolio. We already have the historic return series stored from earlier.

mu = portfolio_returns.mean()
sigma = portfolio_returns.std()

print(f'Our portfolio mean return value is {round(mu*100,2)}%')
print(f'Our portfolio standard deviation value is {round(sigma*100,2)}%')

Our portfolio mean return value is 0.04% 
Our portfolio standard deviation value is 0.85%

Now we generate the necessary draws from the normal distribution with mean 0.04% and standard deviation 0.85%

portfolio_mc = pd.DataFrame([(np.random.normal(loc=mu, scale=sigma, size=252)+1) for x in range(1000)]).T.cumprod()
portfolio_mc.plot(figsize=figsize, legend=False, linewidth=1, alpha=0.2, color='green')

Let us finally now run the Monte Carlo simulation approach but this time create random draws from each individual asset distribution and then construct our portfolio and see if there is any difference in outcome.

for asset in (asset_universe_returns.mean() * 100).round(2).index:
    print(f'The mean return for {asset} is {(asset_universe_returns.mean() * 100).round(2)[asset]}%')
    
print('\n')
for asset in (asset_universe_returns.std() * 100).round(2).index:
    print(f'The mean return for {asset} is {(asset_universe_returns.std() * 100).round(2)[asset]}%')

The mean return for ^DJI is 0.04% 
The mean return for ^IXIC is 0.05% 
The mean return for ^GSPC is 0.04% 
The mean return for ^STOXX50E is 0.01% 
The mean return for ^N225 is 0.04% 
The mean return for ^GDAXI is 0.04% 

The standard deviation for ^DJI is 0.88% 
The standard deviation for ^IXIC is 1.06% 
The standard deviation for ^GSPC is 0.92% 
The standard deviation for ^STOXX50E is 1.25% 
The standard deviation for ^N225 is 1.29% 
The standard deviation for ^GDAXI is 1.19% 

Create our DataFrames of simulated asset returns for each individual asset and store them in a list.

asset_returns_dfs = []
for asset in asset_universe_returns.mean().index:
    mu = asset_universe_returns.mean()[asset]
    sigma = asset_universe_returns.std()[asset]
    asset_mc_rets = pd.DataFrame([(np.random.normal(loc=mu, 
                    scale=sigma, size=252)) for x in range(1000)]).T
    
    asset_returns_dfs.append(asset_mc_rets)

Use list comprehension to iterate through the list of asset return DataFrames and divide the values by the number of assets to represent an equally weighted portfolio.

weighted_asset_returns_dfs = [(returns_df / len(tickers)) for returns_df in asset_returns_dfs]

Add together the DataFrame values using the “reduce” function from the “functools” library (great library by the way, along with the “itertools” library. Loads of great useful functions and definitely worth checking out)

portfolio_constituents_mc = (reduce(lambda x, y: x + y,weighted_asset_returns_dfs) + 1).cumprod()

Finally we plot the resulting Monte Carlo portfolio value simulations.

portfolio_constituents_mc.plot(figsize=figsize, legend=False, linewidth=1, alpha=0.2, color='orange')

We can see immediately that something looks different!! Or well, maybe not immediately but something should jump out at us. If you notice, all the previous simulations, whether Bootstrapping or Monte Carlo have all produce simulations which fall within the ending value bounds of around 0.8 to 1.6. But in the last plot we see these bounds have tightened to around 0.9 to 1.3.

That’s a significant difference, and one that can not be put down to the effects of randomness alone. If you re-run all these simulations a a few times, you will see that the outcomes remain similar, and that the last method will pretty much always produce a tighter range of ending values.

So why is this you may be asking!

Well remember when I mentioned the effects of correlation between individual assets, and the fact that we had to be careful to capture this effect when running our simulations? It is only the last method which fails to capture this effect of correlation.

Let’s take a quick look at the historic correlations between the constituent assets’ returns in the “asset universe” we chose.

NOTE – it is important to calculate the correlation between the asset’s RETURNS, NOT their prices (perhaps more on that in a future article).

Let’s create a nice looking little correlation heatmap to have a look.

ax, fig = plt.subplots(figsize=(12,10))
sns.heatmap(asset_universe_returns.corr(),annot=True)
plt.plot()

Looking at the above values, we can see that all of the assets are positively correlated to some degree, some more than others, bit importantly the values are all positive.

This goes a long way to explaining why our last plot and out last simulation method (the parametric Monte Carlo simulations on constituent assets that were then weighted and summed to represent our portfolio) resulted in a more narrow range of ending values.

The logic is reasonably simple – when two assets are correlated, they tend to move in the same direction at the same time – so if one experiences a rise in value, generally so does the other (with some caveats – again, perhaps more in a future article). This results in portfolios which contain positively correlated assets to experience, on average, more extreme values than a portfolio of totally uncorrelated assets, or indeed that of a portfolio of negatively correlated assets.

That is because if all the constituent assets are highly correlated, they will all tend to move up and down at the same time, causing more volatile swings in value.

So how does this relate to our situation and why is it the cause of the difference in outcomes in the last model?

The first 3 approaches all captured the inherent correlation that exists between the constituent assets, whereas our last approach didn’t.

Approach 1 created our portfolio using real historic daily return values that were actually seen to happen in the markets on the same day – so the moves that were incorporated were real moves that were generated by the underlying process that accounted for and was affected by the real correlation between the assets.

Same logic for approach 2 – although we Bootstrapped returns for our individual assets this time and THEN formed out portfolio, the initial returns that were Bootstrapped were again carefully selected so that all returns in a single draw were taken from the same day for each asset. That again made sure that the values we were extracting were real values that actually took place on a single day, and again were generated by the real underlying process which implicitly accounts for the correlation between the assets.

Now it becomes slightly more subtle as to the difference for approach 3 and 4.

With approach 3, we created our portfolio using real real individual asset returns and THEN ran the parametric Bootstrap process, simulating return series based on the underlying characteristics of the portfolio. Now here is the important part- because the portfolio was FIRST created by using real values of weighted daily returns of individual assets, the price series implicitly accounted for the correlation between the assets. THEN when the Monte Carlo simulation was run, it was being fed input parameters that were calculated on a historic price series that had those correlation relationships implicitly built in. So that method DID capture the effects of correlation.

However, with approach 4 we failed to model the effect of correlation correctly. The individual Monte Carlo simulations for each asset were being fed with parameters that were calculated based on calculated values that were completely independent of each other.

The calculation for the mean and standard deviation of one asset, was being done in a “vacuum” as it were, in a way that was completely independent from the other assets.

With the draws then being made from a normal distribution meant that the individual values drawn on each day for each asset were indeed “random” – that is to to say there was an equal chance that the outcome could be positive or negative for each individual asset, regardless of the output for the other assets.

Well this sounds very much like what we would expect from a series of completely uncorrelated assets – for each to move randomly, irrespective of the moves of the others.

So suddenly we have begun to model the simulated prices series of a basket of uncorrelated assets!!!!

That’s not what we want…so be careful when carrying out these approaches and make sure you are correctly modelling what you are actually trying to model!

Until next time…. 😀