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]:

	Inflation	Nominal Stock Return	Nominal Bond Return
Year
1793	0.072641	-0.076413	-0.072395
1794	0.127533	0.198080	0.153205
1795	0.095134	0.098992	-0.010299
1796	0.006405	-0.035020	-0.060962
1797	-0.035402	0.133346	0.154601
...	...	...	...
2014	-0.000893	0.108304	0.203144
2015	0.013731	-0.047046	-0.059783
2016	0.025000	0.221602	0.057681
2017	0.020705	0.240163	0.097465
2018	0.015512	-0.029244	0.003409

226 rows × 3 columns

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())

Historical average real stock returns: 0.07380060050364495

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__)

yfinance version: 0.2.54

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]:

	CPI	Change
Year
1947	21.480	0.102421
1948	23.680	0.013936
1949	24.010	-0.020825
1950	23.510	0.079541
1951	25.380	0.042159
...	...	...
2021	262.639	0.075781
2022	282.542	0.063403
2023	300.456	0.031079
2024	309.794	0.029994
2025	319.086	NaN

79 rows × 2 columns

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)
)

	Price	Change
Year
2002	34.557056	-0.220025
2003	26.953640	0.379811
2004	37.190929	0.071648
2005	39.855572	0.125383
2006	44.852764	0.141750
2007	51.210640	-0.029341
2008	49.708065	-0.382489
2009	30.695290	0.351860
2010	41.495743	0.242646
2011	51.564503	0.039582
2012	53.605530	0.168366
2013	62.630878	0.225794
2014	76.772575	0.130532
2015	86.793846	-0.027279
2016	84.426178	0.218880
2017	102.905396	0.252306
2018	128.869064	-0.022549
2019	125.963135	0.203136
2020	151.550812	0.207456
2021	182.990982	0.184618
2022	216.774399	-0.083999
2023	198.565521	0.192004
2024	236.690826	0.261561
2025	298.600006	NaN

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]:

	Price_BLV	Change_BLV	Price_VCLT	Change_VCLT	Price	Change
Year
2008	35.901684	0.016191	NaN	NaN	NaN	NaN
2009	36.482979	0.090712	NaN	NaN	NaN	NaN
2010	39.792435	0.063808	38.102699	0.089186	38.947567	0.076497
2011	42.331512	0.257824	41.500919	0.205286	41.916216	0.231555
2012	53.245609	0.031417	50.020470	0.065415	51.633039	0.048416
2013	54.918442	-0.015759	53.292538	0.009313	54.105490	-0.003223
2014	54.052979	0.220951	53.788864	0.186349	53.920921	0.203650
2015	65.996017	-0.078047	63.812355	-0.102151	64.904186	-0.090099
2016	60.845230	0.046441	57.293850	0.113832	59.069540	0.080137
2017	63.670952	0.077371	63.815746	0.103114	63.743349	0.090243
2018	68.597244	-0.002280	70.396049	-0.022199	69.496647	-0.012240
2019	68.440842	0.230346	68.833321	0.236660	68.637081	0.233503
2020	84.205948	0.067297	85.123421	0.060251	84.664684	0.063774
2021	89.872719	-0.042398	90.252220	-0.039607	90.062469	-0.041003
2022	86.062271	-0.181506	86.677612	-0.155396	86.369942	-0.168451
2023	70.441437	-0.012861	73.208260	0.024118	71.824848	0.005629
2024	69.535522	-0.018190	74.973869	-0.006960	72.254696	-0.012575
2025	68.270660	NaN	74.452049	NaN	71.361355	NaN

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]:

	Real_Stock_Return	Real_Bond_Return
Year
1793	-0.138960	-0.135214
1794	0.062568	0.022768
1795	0.003523	-0.096274
1796	-0.041162	-0.066938
1797	0.174941	0.196976
...	...	...
2014	0.109295	0.204220
2015	-0.059954	-0.072518
2016	0.191806	0.031883
2017	0.215006	0.075203
2018	-0.044072	-0.011919

226 rows × 2 columns

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

Last year in McQuarrie data: 2018

Out[17]:

	Inflation	Nominal_Stock_Returns	Nominal_Bond_Returns	Real_Stock_Return	Real_Bond_Return
Year
2019	0.025998	0.203136	0.233503	0.172650	0.202248
2020	0.013553	0.207456	0.063774	0.191310	0.049549
2021	0.075781	0.184618	-0.041003	0.101170	-0.108557
2022	0.063403	-0.083999	-0.168451	-0.138614	-0.218030
2023	0.031079	0.192004	0.005629	0.156074	-0.024684
2024	0.029994	0.261561	-0.012575	0.224824	-0.041330
2025	NaN	NaN	NaN	NaN	NaN

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.