Backtesting the US market with ca. 100 years of data¶
This page by Damodaran, with curated data since 1928, motivated me to combine all possible portfolios and get statistics.
The 6 asset classes are:
- S&P 500 including dividends
- Cash (i.e., short-terms US Treasury bills)
- Long-term US government bonds (10y)
- Invest Grade Corporate Bonds (Baa)
- Real Estate
- Gold
In [1]:
import numpy as np
import pandas as pd # Note: needs xlrd and openpyxl to be installed to read excel files
from tqdm import tqdm
import datetime
import ssl
import matplotlib.pyplot as plt
import plotly.express as px
import plotly.graph_objects as go
import plotly.io as pio
pd.set_option('plotting.backend', 'plotly')
pio.renderers.default = 'notebook'
ssl._create_default_https_context = ssl._create_unverified_context # Needed to download data
In [2]:
def cumulative_return(x):
"""Cumulative return to be used with df.rolling().apply().
Consider that x is a one-period return, for example:
[ 0.5, 0.2, -0.1] means +50%, +20%, -10%, and what this does is total_return = (1 + 0.5) * (1 + 0.2) * (1 - 0.1) - 1
"""
return (x + 1).prod() - 1
def drawdown(x):
if isinstance(x, np.ndarray):
x = pd.Series(x)
cumulative = (1 + x).cumprod() # Cumulative returns
peak = np.maximum.accumulate(cumulative) # equivalent of cumulative.cummax() which is only available in pandas and not numpy
return (cumulative / peak) - 1 # Drawdown
def max_drawdown(x):
"""Max drawdown to be used with df.rolling().apply().
This is the largest loss that an investor that look at his portfolio at each point in time would have seen.
"""
return drawdown(x).min()
# Let's see some examples of max drawdown
for returns in (
[0.1, 0.2, 0.3], # expected: 0
[0.1, 0.2, 0.3, -0.1], # expected: -0.1
[0.2, 0.3, -0.4, -0.3, 0.2], # expected: 1 - (1-0.4)*(1-0.3) = 1 - 0.42 = -0.58, NOTE that this is the drop observed from the highest value of the asset (1.2*1.3 = 1.56 of the starting value) considered as unity.
# This is different from the drop since the start: cumulative_return([0.2, 0.3, -0.4, -0.3]) = -0.34
[0.1, -0.05, 0.2, -0.1, 0.05, 0.3, -0.2, 0.1, -0.15, 0.05], # there are several drawdowns, -0.05, -0.1, -0.2, 1 - (1-0.2)*1.1*(1-0.15) = -0.25 where in the last the value had not enough bump up to recover the previous losses
):
plt.figure(figsize=(10,1))
plt.plot(drawdown(pd.Series(returns)), marker='o')
plt.xticks(range(len(returns)), range(len(returns)))
plt.xlabel(f"Max Drawdown: {max_drawdown(pd.Series(returns)):.2f}")
plt.show()
In [3]:
# Download the data
url_dataset = "https://www.stern.nyu.edu/~adamodar/pc/datasets/histretSP.xls"
SKIP_ROWS = 19
START_YEAR = 1928
LAST_YEAR = datetime.date.today().year -1
dfy = (
pd.read_excel(url_dataset, sheet_name="Returns by year", skiprows=SKIP_ROWS, nrows=LAST_YEAR-START_YEAR+1)
.iloc[:, :7] # keep only the first 7 columns
.set_index("Year")
.rename(columns={
'S&P 500 (includes dividends)': 'S&P500',
'3-month T.Bill': "Cash",
'US T. Bond (10-year)': "TBond",
' Baa Corporate Bond': "CorpBond",
'Real Estate': "RealEstate",
'Gold*': "Gold"
})
)
dfy
Out[3]:
In [4]:
fig = px.line(dfy*100, title="Annual returns", labels={"value": "Return (%)", "Year": "Year"})
fig.update_layout(title="Timeline of 1y returns", yaxis_title="Return (%)", margin=dict(l=5, r=5, t=40, b=5), width=1000)
fig.show()
fig = go.Figure()
for col in dfy.columns:
fig.add_trace(go.Violin(y=dfy[col]*100, name=col))
fig.update_layout(title="Distribution of 1y returns", yaxis_title="Return (%)", margin=dict(l=5, r=5, t=40, b=5), width=1000)
fig.show()
In [5]:
years = [3, 5, 10, 20, 30]
for year in years:
df = dfy.copy() # never modify the original data in dfy
fig = go.Figure()
for asset in df.columns:
col = f'{asset}-{year}y'
df[col] = df[asset].rolling(year).apply(cumulative_return, raw=True)
fig.add_trace(go.Violin(y=df[col]*100, name=asset))
fig.update_layout(title=f"Distribution of RETURNS in a {year}y window", yaxis_title="Return (%)", margin=dict(l=5, r=5, t=40, b=5), width=1000, height=200)
fig.show()
In [6]:
# Same with the drawdowns (NOTE: data are yearly, so the drawdown in the minimum return of the time window at the end of the year)
years = [3, 5, 10, 20, 30]
for year in years:
df = dfy.copy() # never modify the original data in dfy
fig = go.Figure()
for asset in df.columns:
col = f'{asset}-{year}y'
df[col] = df[asset].rolling(year).apply(max_drawdown, raw=True)
fig.add_trace(go.Violin(y=df[col]*100, name=asset))
fig.update_layout(title=f"Distribution of MAX-DRAWDOWNS in a {year}y window", yaxis_title="Return (%)", margin=dict(l=5, r=5, t=40, b=5), width=1000, height=200)
fig.show()
In [7]:
years_of_custody = 20
custom_portfolio ={
"S&P500": 0.2,
"Cash": 0.1,
"TBond": 0.2,
"CorpBond": 0.1,
"RealEstate": 0.2,
"Gold": 0.2
}
assert sum(custom_portfolio.values()) == 1, "The sum of the weights must be 1"
df = dfy.copy() # never modify the original data in dfy
# Compute the average return and the max drawdown for the custom portfolio
df["Custom"] = 0
for asset, weight in custom_portfolio.items():
df["Custom"] += df[asset] * weight
average_return = df["Custom"].mean() # it is a mean, it is irrelevant the years of custody
df["Custom-RET"] = df["Custom"].rolling(years_of_custody).apply(cumulative_return, raw=True)
annualized_min_return = (1 + df["Custom-RET"].min())**(1/years_of_custody) - 1
annualized_max_return = (1 + df["Custom-RET"].max())**(1/years_of_custody) - 1
print(f"Return min annualized: {annualized_min_return:.2%}")
print(f"Return average: {average_return:.2%}")
print(f"Return max annualized: {annualized_max_return:.2%}")
df["Custom-MDD"] = df["Custom"].rolling(years_of_custody).apply(max_drawdown, raw=True)
print(f"Max Drawdown: {df['Custom-MDD'].min():.2%}")
Make a random generator of portfolios¶
Otherwise in grid search with 6 assets and 10% step I would have 1e6 combinations, 5% 64e6 combinations.
In [ ]:
In [8]:
def generate_portfolios(assets, n_portfolios, step=0.1):
portfolios_list = []
# Add for each asset 100% one asset (would be inefficient to have it in the loop)
# (also keeping the order of columns as in the list of assets)
for asset in assets:
portfolios_list.append({asset: 1.0})
# My cash portfolio
for i in range(1, int(1/step)):
portfolios_list.append({"Cash": 1-step*i, "S&P500": step*i})
for _ in range(10): # Stop if you can do all n_portfolios I asked in 10 tries
for _ in range(n_portfolios):
np.random.shuffle(assets) # Shuffle assets to randomize the order
portfolio_weights = {}
total_weight = 0.0
for asset in assets:
if total_weight < 1.0:
max_possible_weight = 1.0 - step - total_weight
steps_possible = int(max_possible_weight * 10)
weight = np.random.randint(0, steps_possible + 1) / 10
portfolio_weights[asset] = weight
total_weight += weight
# Adjust the last weight if the total is not exactly 1.0 due to rounding
if not np.isclose(total_weight, 1.0):
last_asset = assets[-1]
portfolio_weights[last_asset] += 1.0 - total_weight
portfolios_list.append(portfolio_weights)
# Convert the list to DataFrame and ensure unique portfolios
portfolios_df = pd.DataFrame(portfolios_list).drop_duplicates(keep="first").fillna(0)
# If the DataFrame has more portfolios than needed, trim the excess
if len(portfolios_df) > n_portfolios:
portfolios_df = portfolios_df.iloc[:n_portfolios]
break
# if the DataFrame has less portfolios than needed, try again but get rid of duplicates
portfolios_list = portfolios_df.to_dict(orient='records')
return portfolios_df.reset_index(drop=True)
In [9]:
df_portfolios = generate_portfolios(list(custom_portfolio.keys()), 500)
df_portfolios.head(20)
Out[9]:
In [10]:
df_portfolios.sum(axis=1).value_counts()
Out[10]:
In [11]:
# give me the distbution histogram of each column
fig = go.Figure()
for col in df_portfolios.columns:
fig.add_trace(go.Histogram(x=df_portfolios[col], name=col))
fig.update_layout(title="Distribution of weights in random portfolios", xaxis_title="Weight", yaxis_title="Count", margin=dict(l=5, r=5, t=40, b=5), width=1000)
fig.show()
In [12]:
years_of_custody = 10
assets = list(custom_portfolio.keys())
df = dfy.copy()
df_portfolios["Return"] = 0.0
df_portfolios["MaxDrawdown"] = 0.0
for i, weights in tqdm(enumerate(df_portfolios[assets].to_numpy()), total=len(df_portfolios)):
portfolio_return = (df[assets] * weights).sum(axis=1) # pd.Series
df_portfolios.at[i, "Return"] = portfolio_return.mean()
df_portfolios.at[i, "MaxDrawdown"] = portfolio_return.rolling(years_of_custody).apply(max_drawdown, raw=True).min()
df_portfolios
Out[12]:
In [13]:
df_portfolios.head(20)
Out[13]:
In [15]:
# scatter plot Returns vs MaxDrawdown, including all columns in the hover
fig = px.scatter(df_portfolios.iloc[::-1], x="MaxDrawdown", y="Return", hover_data=df_portfolios.columns, title="Random portfolios")
fig.update_layout(title="Random portfolios", margin=dict(l=5, r=5, t=40, b=5), width=1000)
fig.update_yaxes(tickformat=".0%")
fig.update_xaxes(tickformat=".0%")
fig.data[0].marker.color = ['blue'] * (len(df_portfolios) - 6 -9) + ['green'] * 9 + ['red'] * 6 # Adjust 'blue' to your default color
fig.data[0].marker.size = [7] * (len(df_portfolios) - 9 - 6) + [15]*(9+6) # Adjust 5 to your default size
fig.show()
Conclusions¶
I'm not happy with this: I'm not sure if the metrics I'm using (mean returns, max drawdown) are the best since are not much dependent on the custody-period-length.
It is clear that an optimized portfolio does better than a simple cash+s&p adjusted for risk-tolerance, but what is the proper way to condider the fact that my life-experience is not ergodic and I don't want any chance of tail risk?