Bond and Stocks lead to similar results in the long run¶
Data are from the US stock market, since 1793: 250y of data.
The thesis is that when you take long term bonds, except for some periods (1940-1980) bonds and stocks surprisingly lead to similar results!
- Source https://www.edwardfmcquarrie.com/?p=579
- Paper: https://dx.doi.org/10.2139/ssrn.3805927
- backup @ my Drive/Stuff
- Bonds are intended as 15+ years maturity (long term) and a mix of government and corporate bonds.
In [1]:
import numpy as np
import pandas as pd # Note: needs xlrd and openpyxl to be installed to read excel files
import ssl
import plotly.express as px
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]:
import pandas as pd
import requests
import io
url_dataset = "http://www.edwardfmcquarrie.com/wp-content/uploads/2021/07/Real-returns-on-stocks-and-bonds-1793-to-2019-version-2-0.xlsx"
headers = {
'User-Agent': 'Mozilla/5.0 (Macintosh; Intel Mac OS X 10_15_7) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/91.0.4472.114 Safari/537.36'
}
response = requests.get(url_dataset, headers=headers)
response.raise_for_status() # Raise an HTTPError for bad responses
# Read the content into a DataFrame
dfy = (
pd.read_excel(io.BytesIO(response.content), sheet_name="real returns 1793-2019", skiprows=None, nrows=None)
.dropna(subset="Nominal stock return") # remove first empty rows
.rename(columns={
'To January of:': 'Year',
'Annual inflation relative': 'Inflation',
'Nominal stock return': 'Nominal Stock Return',
'Nominal bond return': 'Nominal Bond Return',
})
[['Year', 'Inflation', 'Nominal Stock Return', 'Nominal Bond Return']]
.assign(Year=lambda x: x["Year"] - 1) # refer to the retunr of the year, not at Jenuary of the next year
.assign(Inflation= lambda x: x["Inflation"] -1 ) # get inflation as pct change
.set_index("Year")
)
dfy
Out[2]:
In [3]:
dfy.plot().update_layout(yaxis_title="Percentage Change", yaxis_tickformat=".0%" ,margin=dict(l=5, r=5, t=40, b=5), width=1000)
In [4]:
# for each decade calculate the average return
# add years between 1790 and 1793
dfy_decade = dfy.copy()
dfy_decade.loc[1790:1792] = np.nan
dfy_decade = dfy_decade.groupby(dfy.index // 10).mean()
dfy_decade.index = dfy_decade.index * 10
px.bar(dfy_decade, barmode='group').update_layout(
title="Average Returns by Decade",
yaxis_title="Percentage Change", yaxis_tickformat=".0%" ,margin=dict(l=5, r=5, t=40, b=5), width=1000)
In [5]:
dfy_cum = dfy + 1
dfy_cum = dfy_cum.cumprod()
dfy_cum.plot().update_layout(yaxis_type="log").update_layout(
title="Nominal Cumulative Returns",
yaxis_title="Value relative to first year",margin=dict(l=5, r=5, t=40, b=5), width=1000
).add_vrect(x0=1940, x1=1980, fillcolor="green", opacity=0.1, layer="below", line_width=0
)
In [6]:
# Real returns
dfy_real = (
dfy
.assign(Real_Stock_Return=lambda x: (1 + x["Nominal Stock Return"]) / (1 + x["Inflation"]) - 1)
.assign(Real_Bond_Return=lambda x: (1 + x["Nominal Bond Return"]) / (1 + x["Inflation"]) - 1)
[['Real_Stock_Return', 'Real_Bond_Return']]
)
dfy_real_cum = dfy_real + 1
dfy_real_cum = dfy_real_cum.cumprod()
dfy_real_cum.plot().update_layout(yaxis_type="log").update_layout(
title="Real Returns",
yaxis_title="Value relative to first year",margin=dict(l=5, r=5, t=40, b=5), width=1000
).add_vrect(x0=1940, x1=1980, fillcolor="green", opacity=0.1, layer="below", line_width=0
)
In [7]:
print("Historical average real stock returns:", dfy_real.Real_Stock_Return.mean())
In [8]:
# let's alling the 1940-1980 period
dfy_real_mod = dfy_real.copy()
dfy_real_mod.loc[1940:1980, "Real_Bond_Return"] = 0.0738
dfy_real_cum = dfy_real_mod + 1
dfy_real_cum = dfy_real_cum.cumprod()
dfy_real_cum.plot().update_layout(yaxis_type="log").update_layout(
title="Real Returns (1940-1980 Bond Returns MODIFIED from flat to 7.4%)",
yaxis_title="Value relative to first year",margin=dict(l=5, r=5, t=40, b=5), width=1000
).add_vrect(x0=1940, x1=1980, fillcolor="green", opacity=0.1, layer="below", line_width=0
)
Conclusions¶
The thesis of McQuarrie is that one can say that stocks' returns outperform (long term) bonds' returns (Jeremy Siegel's thesis) only if you look at data from the 20th century and later.
Considering 250 years of US data, the returns of the two assets are surprisingly similar
the only multi-decade exception was the 1940-1980 period, where real return on bonds stayed flat (while the real return on stocks was about the same as the historical average).
Conceptually, this should not be very surprising as equity and long-term debt can be seen as the same risk of investment, even if it is a common conception that stocks should lead to higher returns as they can be wiped out sooner than bonds.
Follow-up¶
- Read McQuarrie's paper, to understand if my conceptual explanation is the same of of his
- read his explanation on why the 1940-1980 period was an exception.
Part 2: fetch old data with new data that can be updated¶
From the Appendix of McQuarrie's paper, I can see that he used data from the following sources:
- Stocks: 1926 to present - CRSP total market index (NYSE only until 1962)
- Bonds: 1974 to present - Uses the SBBI long corporate bond index (20y) rather than the SBBI long government index used by Siegel. This adds just under 20 basis points to the annualized returns. Rationale: these top-grade bonds (Aaa, Aa) provide a proxy for the entire investment grade space, if their returns are taken as near the midpoint of Treasury, agency, mortgage, and medium grade corporate bond returns (A, Baa).
- Inflation: 1913 to present - BLS CPI-U index (CPI before 1978). Uses January values of the CPI-U downloaded from the BLS website.
I don't have these sources at my disposal, and I want to use something more practical like EFTs and FRED data.
In [9]:
import os
import dotenv # pip install python-dotenv
import pandas as pd
import yfinance as yf # pip install yfinance
from fredapi import Fred # pip install fredapi
dotenv.load_dotenv()
print("yfinance version:", yf.__version__)
In [10]:
fred = Fred(api_key=os.getenv("FRED_API_KEY"))
# Fetch historical data on US inflation (CPI)
cpi_data = fred.get_series('CPIAUCSL')
cpi_df = pd.DataFrame(cpi_data, columns=['CPI'])
cpi_df.index.name = 'Date'
cpi_year_df = (
cpi_df
.loc[cpi_df.index.month == 1]
.assign(
Year=lambda df: df.index.year,
Change=lambda df: df["CPI"].pct_change().shift(-1)
)
.set_index("Year")
)
cpi_year_df
Out[10]:
In [11]:
# compare inplotly express cpi_year_df.Change vs dfy.Inflation
df_infl = pd.merge(
left=cpi_year_df.Change,
right=dfy.Inflation,
left_index=True,
right_index=True,
how="outer"
).rename(columns={"Change": "From FRED", "Inflation": "From McQuarrie"})
df_infl.plot(
).update_traces(
selector=dict(name="From FRED"), line=dict(width=5)
).update_layout(
title="Inflation Rates",
yaxis_title="Percentage Change", yaxis_tickformat=".0%",
margin=dict(l=5, r=5, t=40, b=5), width=1000
)
In [12]:
yfstocks_df = (
yf.Ticker("VTI") # Vanguard Total Stock Market ETF (US stocks)
.history(period="max", interval="1mo", auto_adjust=True) # adjust=True to account for splits and dividends
[["Close"]]
.rename(columns={"Close": "Price"})
.query("index.dt.month == 1")
.assign(
Year=lambda df: df.index.year,
Change=lambda df: df["Price"].pct_change().shift(-1)
)
.set_index("Year")
)
display(yfstocks_df)
df_stocks = pd.merge(
left=yfstocks_df.Change,
right=dfy["Nominal Stock Return"],
left_index=True,
right_index=True,
how="outer"
).rename(columns={"Change": "From YF", "Nominal Stock Return": "From McQuarrie"})
px.scatter(
df_stocks.reset_index(),
x="From YF",
y="From McQuarrie",
hover_data=["Year"]
).update_layout(
title="Stock Returns: compare MCQuarrie-Nominal vs VTI (YahooFin)",
xaxis_title="Percentage Change (Yahoo Finance VTI)",
yaxis_title="Percentage Change (McQuarrie)",
margin=dict(l=5, r=5, t=40, b=5), width=600, height=500
).add_shape(
type='line',
x0=-0.4, y0=-0.4,
x1=0.4, y1=0.4,
line=dict(color='black', width=1)
)
In [13]:
# Just for curiosity: let's compare with SP500 where I have more data (back to 1985 instead of 2002 for VTI), but it is a different sampling
yfstocks500_df = (
yf.Ticker("^GSPC") # S&P 500 Index
.history(period="max", interval="1mo", auto_adjust=True)
[["Close"]]
.rename(columns={"Close": "Price"})
.query("index.dt.month == 1")
.assign(
Year=lambda df: df.index.year,
Change=lambda df: df["Price"].pct_change().shift(-1)
)
.set_index("Year")
)
df_stocks500 = pd.merge(
left=yfstocks500_df.Change,
right=dfy["Nominal Stock Return"],
left_index=True,
right_index=True,
how="outer"
).rename(columns={"Change": "From YF", "Nominal Stock Return": "From McQuarrie"})
px.scatter(
df_stocks500.reset_index(),
x="From YF",
y="From McQuarrie",
hover_data=["Year"]
).update_layout(
title="Stock Returns: compare MCQuarrie-Nominal vs S&P500 (YahooFin)",
xaxis_title="Percentage Change (Yahoo Finance VTI)",
yaxis_title="Percentage Change (McQuarrie)",
margin=dict(l=5, r=5, t=40, b=5), width=600, height=500
).add_shape(
type='line',
x0=-0.4, y0=-0.4,
x1=0.4, y1=0.4,
line=dict(color='black', width=1)
)
In previous cells I've considered both SP500 (which dates back to 1985) and Vanguard Total Bond Market Index Fund (which dates back to 2002).
It would have been nicer to use the SP500 which has more data but the VTI better matches the data used by McQuarrie:
the returns from these two sets of data are very very similar, so I can assume they are the same.
In [14]:
bond_etfs = [
"BLV", # Vanguard Long-Term Bond ETF (US government bonds)
"VCLT", # Vanguard Long-Term Corporate Bond ETF (US corporate bonds)
]
BLEND_RATIO = 0.5
yfbonds_df = (
yf.Ticker("BLV") # Vanguard Long-Term Bond ETF (US government bonds)
.history(period="max", interval="1mo", auto_adjust=True)
[["Close"]]
.rename(columns={"Close": "Price"})
.query("index.dt.month == 1")
.assign(
Year=lambda df: df.index.year,
Change=lambda df: df["Price"].pct_change().shift(-1)
)
.set_index("Year")
).join(
other=(
yf.Ticker("VCLT") # Vanguard Long-Term Bond ETF (US government bonds)
.history(period="max", interval="1mo", auto_adjust=True)
[["Close"]]
.rename(columns={"Close": "Price"})
.query("index.dt.month == 1")
.assign(
Year=lambda df: df.index.year,
Change=lambda df: df["Price"].pct_change().shift(-1)
)
.set_index("Year")
),
how="outer",
lsuffix="_" + bond_etfs[0],
rsuffix="_" + bond_etfs[1]
).assign(
Price=lambda df: BLEND_RATIO*df["Price_" + bond_etfs[0]] + (1-BLEND_RATIO)*df["Price_" + bond_etfs[1]],
Change=lambda df: BLEND_RATIO*df["Change_" + bond_etfs[0]] + (1-BLEND_RATIO)*df["Change_" + bond_etfs[1]]
)
yfbonds_df
Out[14]:
In [15]:
df_bonds = pd.merge(
left=yfbonds_df.Change,
right=dfy["Nominal Bond Return"],
left_index=True,
right_index=True,
how="outer"
).rename(columns={"Change": "From YF", "Nominal Bond Return": "From McQuarrie"})
px.scatter(
df_bonds.reset_index(),
x="From YF",
y="From McQuarrie",
hover_data=["Year"]
).update_layout(
title="Bonds Returns: compare MCQuarrie-Nominal vs MixedBonds (YahooFin)",
xaxis_title="Percentage Change (Yahoo Finance Bond)",
yaxis_title="Percentage Change (McQuarrie)",
margin=dict(l=5, r=5, t=40, b=5), width=600, height=500
).add_shape(
type='line',
x0=-0.4, y0=-0.4,
x1=0.4, y1=0.4,
line=dict(color='black', width=1)
)
For bonds I obtained a very good agreement by mixing 50:50 the returns of a long term treasury ETF with invetment grade corporate bonds ETF.
TODO: I could optimize further the ratio to improve the agreement but there is not much margin for improvement.
Now, I can update McQuarrie's data with the new data I have obtained.
In [16]:
dfy_real
Out[16]:
In [17]:
last_year_mcquarrie = dfy_real.index.max()
print("Last year in McQuarrie data:", last_year_mcquarrie)
dfy_real_new = (
df_infl[["From FRED"]].rename(columns={"From FRED": "Inflation"}
).join(df_stocks[["From YF"]].rename(columns={"From YF": "Nominal_Stock_Returns"})
).join(df_bonds[["From YF"]].rename(columns={"From YF": "Nominal_Bond_Returns"})
)
.query("index > @last_year_mcquarrie")
.assign(
Real_Stock_Return=lambda x: (1 + x["Nominal_Stock_Returns"]) / (1 + x["Inflation"]) - 1,
Real_Bond_Return=lambda x: (1 + x["Nominal_Bond_Returns"]) / (1 + x["Inflation"]) - 1
)
)
dfy_real_new
Out[17]:
In [18]:
# let's alling the 1940-1980 period
dfy_real_mod = dfy_real.copy()
dfy_real_mod.loc[1940:1980, "Real_Bond_Return"] = 0.0738
dfy_real_mod = pd.concat([dfy_real_mod, dfy_real_new[['Real_Stock_Return', 'Real_Bond_Return']]])
dfy_real_cum = dfy_real_mod + 1
dfy_real_cum = dfy_real_cum.cumprod()
dfy_real_cum.plot().update_layout(yaxis_type="log").update_layout(
title=f"Real Returns (1940-1980 Bond Returns MODIFIED from flat to 7.4%), w/post-{last_year_mcquarrie} data",
yaxis_title="Value relative to first year",margin=dict(l=5, r=5, t=40, b=5), width=1000
).add_vrect(
x0=1940, x1=1980, fillcolor="green", opacity=0.1, layer="below", line_width=0
).add_vrect(
x0=2017, x1=2024, fillcolor="purple", opacity=0.1, layer="below", line_width=0
)
Surprising: there is a new divergence between stocks (outperforming) and bonds (underperforming)!
Is this a new 1940-moment where bonds will underperform for decades?
DISCLAIMER: re-check all calculations and data.