Leverage ETFs¶
Goal is to calculate historical returns for a leveraged SP500 ETF, since there is no historical data for e.g. the UPRO
We try to hit the UPRO (3x SP500) as a benchmark and use FED Fund Rate and SOFR as proxy for borrowing cost.
Step 1: Derive SP500 returns with dividends included (benchmark: SP500 Total Return).
Step 2: Calculate leveraged returns (benchmark: UPRO).
Step 3: Analyze return and drawdown for single and reoccuring investment.Results:
Cost of UPRO is 1.19% p.a. higher than expected
Drawdown of ~90% realized on long holding periods
1. Import Data¶
In [749]:
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
def read_fed(path, date_col = "observation_date", value_col = "value"):
df = pd.read_csv(path, index_col=None)
df[date_col] = pd.to_datetime(df[date_col], format = "%Y-%m-%d")
value_col_old = df.columns[1]
df.rename(columns={value_col_old: "value"}, inplace=True)
return df
def read_yfinance(path, date_col = "Date", value_col = "Close"):
df = pd.read_csv(path, index_col=None)
df[date_col] = pd.to_datetime(df[date_col].astype(str).str.split(" ").str[0], format = "%Y-%m-%d")
return df[[date_col, value_col]]
# treasury for maturities 1, 3, 5, 7, 10, 20, 30 years
treasury = {}
maturities = [1, 3, 5, 7, 10, 20, 30]
for maturity in maturities:
treasury[maturity] = read_fed(f"data/DGS{maturity}.csv", date_col = "observation_date")
# read S&P 500 price
sp500 = read_yfinance("data/SP500.csv", date_col = "Date", value_col = "Close")
# read SP500 dividend yield (percentage of 12 month dividend per share)
sp500_div = pd.read_csv("data/SP500_dividends.csv", sep = ";", decimal = ".", index_col=None)
sp500_div["Date"] = pd.to_datetime(sp500_div["Date"], format = "%b %d, %Y")
sp500_div["Value"] = pd.to_numeric(sp500_div["Value"].str.replace("%","").str.strip(), errors='coerce')
# read SP500 total return index (for validation)
sp500_tr = read_yfinance("data/SP500_TR.csv", date_col = "Date", value_col = "Close")
# read UPRO (for validation)
upro = read_yfinance("data/UPRO.csv", date_col = "Date", value_col = "Close")
In [750]:
plt.figure()
sp500.set_index("Date")["Close"].apply(np.log).plot(
title="S&P 500 (log) Price",
ylabel="(log) Price",
xlabel="Date",
)
plt.figure()
sp500_div.set_index("Date")["Value"].plot(
title="S&P 500 Dividend Yield (12 month dividend per share)",
ylabel="Dividend Yield (%)",
xlabel="Date",
)
Out[750]:
<Axes: title={'center': 'S&P 500 Dividend Yield (12 month dividend per share)'}, xlabel='Date', ylabel='Dividend Yield (%)'>
2. Calculate Leveraged Daily Returns¶
2.1 Calculate Daily Returns $r_{basis}(t)$¶
To easen the calculation, we fill every date for which no price is given with the last available price.
In practice dividens are paid as a lumb sum at the end of each month. In order to calculate daily returns, we assume that the monthly dividend $y_m$ is paid daily and is equally distributed.
$$ y(t) = \left(1 + y_m\right)^{\frac{1}{365}} - 1 \text{ for all } t\in m$$
Daily returns are sum of daily price returns and daily dividend returns. We test the results against the SP500 Total Return Benchmark.The residuals inbetween the 1th and 99th quantile are low. The density of return deltas forms a left tail, since benchmark books a whole chunk of dividends at once.
$$ r_{basis}(t) = \frac{P(t)}{P(t-1)}-1 + \frac{y(t)}{365}$$
In [751]:
# get prices and dividends
prices = sp500.set_index("Date").rename(columns={"Close": "og_price"})
dividends_og = sp500_div.set_index("Date").rename(columns={"Value": "og_div_yield"})
full_date_range = pd.date_range(sp500["Date"].min(), sp500["Date"].max(), freq="D")
prices_and_dividends = prices.merge(dividends_og, left_index=True, right_index=True, how="outer").reindex(full_date_range).sort_index()
# calculate daily return including dividend yield
prices_and_dividends["price_intpld"] = prices_and_dividends["og_price"].ffill()
prices_and_dividends["div_yield_intpld"] = prices_and_dividends["og_div_yield"].ffill() / 100
prices_and_dividends["div_yield_daily"] = (1 + prices_and_dividends["div_yield_intpld"])**(1/365) - 1
prices_and_dividends["return_daily"] = (1 + prices_and_dividends["price_intpld"].pct_change()) * (1+ prices_and_dividends["div_yield_daily"])
# set starting value to original price and accumulate returns
prices_and_dividends.iloc[0, prices_and_dividends.columns.get_loc("return_daily")] = prices_and_dividends["og_price"].iloc[0]
prices_and_dividends["SP500"] = prices_and_dividends["return_daily"].cumprod()
# output as csv
prices_and_dividends.to_csv(r"output\SP500_prices_and_dividends.csv")
# plot SP500 Price vs SP500 price with dividends
prices_and_dividends[["SP500", "og_price"]].rename(columns={"og_price": "SP500", "SP500": "SP500 incl. div."}).apply(np.log).plot(
title="SP500 incl. dividends vs SP500",
ylabel="Cumulative Return (log)",
xlabel="Date"
)
## validate against benchmark SP500_TR
# merge with S&P 500 total return index
benchmark_dates = pd.date_range(sp500_tr["Date"].min(), sp500_tr["Date"].max(), freq="D")
benchmark = sp500_tr.set_index("Date").reindex(full_date_range).ffill().rename(columns={"Close": "og_SP500_TR"})
benchmark["benchmark_daily_return"] = 1 + benchmark.pct_change()
# rebase SP500_TR to match starting value of SP500 incl. div
benchmark.loc[sp500_tr["Date"].min(), "benchmark_daily_return"] = prices_and_dividends.loc[sp500_tr["Date"].min(), "SP500"]
benchmark["SP500_TR"] = benchmark["benchmark_daily_return"].cumprod()
prices_and_dividends = prices_and_dividends.merge(benchmark[["benchmark_daily_return", "SP500_TR"]], left_index=True, right_index=True)
# calc delta returns
plt.figure()
prices_and_dividends.loc[benchmark_dates, ["SP500", "SP500_TR"]].rename(columns={"SP500": "SP500 incl. div"}).apply(np.log).plot(
title="S&P 500 Estimate vs. Total Return Benchmark",
ylabel="Level (log)",
xlabel="Date",
linewidth=0.8
)
# plot return deltas
delta_returns = prices_and_dividends["return_daily"] - prices_and_dividends["benchmark_daily_return"]
q99 = delta_returns.quantile(0.99)
q01 = delta_returns.quantile(0.01)
plt.figure()
filtered_data = delta_returns[(delta_returns <= q99) & (delta_returns >= q01)]
filtered_mean = filtered_data.mean()
ax = filtered_data.plot(
title="SP500 incl. dividends returns - SP500_TR returns (1st–99th Percentiles)",
ylabel="Daily return difference (decimal)",
xlabel="Date"
)
ax.axhline(filtered_mean, color = "red", linestyle = "--", label = f"Mean = {filtered_mean:.6f}")
ax.legend()
# plot return delta density
plt.figure()
filtered_data.plot(
kind = "density",
title="Density of Daily Return Differences (1st–99th Percentiles)",
xlabel="Daily return difference (decimal)",
)
Out[751]:
<Axes: title={'center': 'Density of Daily Return Differences (1st–99th Percentiles)'}, ylabel='Density'>
<Figure size 640x480 with 0 Axes>
2.2 Calculate Leveraged daily returns¶
Sounds easy in theory. Calculate daily profits of underlying and multiply with leverage. Subtract cost for the ETF (Expense Ratio ER(t)) and cost of borrowing money for the leverage (Borrowing Cost BC(t)). Also subtract a daily Adjustment Erorr A(t). We will set the adjustment error to 3.09% p.a. for all t as calculated in Section 2.5
$$r_{LETF}(t) = L \cdot r_{basis}(t) - ER(t) - (L-1)\cdot BC(t) - A(t)$$
Cost of borrowing money is given by overnight borrowing rate. We use FED interest rate starting from 1955 as proxy for overnight borrowing rate until SOFR is introduced in 2018. For 1932 to 195 we use the Low Value of the Federal Funds Rate for the Indicated Date Published in The Wall Street Journal, even though this dataset is messy.
$$ BC(t) = \text{SOFR}(t)\cdot \mathbb{1}_{t > 2018} + \text{FED IR}(t)\cdot \mathbb{1}_{t <= 2018}$$
The total expense ratio for the UPRO is around 0.89% p.a. (regarding this factsheet). In reality UPRO also pays a swap spread. To approximate conservatively, we choose ER to be 1%.
$$ ER(t) = (1.001)^\frac{1}{365} - 1$$
The Adjustment factor is set to 2$ p.a. to fit the curve of the UPRO better. It seems that after 2024, the cost of the UPRO increased. adjustment-factor + expense-ratio add up to 3$ p.a. which is much higher than the anticipated 0.91%.
We calculate daily leveraged returns for all $L$ in $\{1, 1.25, ..., 3.75, 4\}$ and compare the 3x leverage with the UPRO.
In [752]:
returns = pd.read_csv(r"output\SP500_prices_and_dividends.csv", index_col = 0)
returns["return_daily"] = returns["return_daily"] - 1
returns.index = pd.to_datetime(returns.index)
# federal funds rate (monthly)
fed_funds = read_fed("data/FEDFUNDS.csv", date_col = "observation_date").set_index("observation_date").rename(columns={
"value": "FED eff. rate"
})
# SOFR rate
sofr = read_fed("data/SOFR.csv", date_col = "observation_date").set_index("observation_date").rename(columns={
"value": "SOFR"
})
# Low Value of the Federal Funds Rate for the Indicated Date Published in The Wall Street Journal for before 1955
ffwsjlow = read_fed("data/FFWSJLOW.csv", date_col = "observation_date").set_index("observation_date").rename(columns={
"value": "ffwsjlow"
})
# daily index over the whole period
full_date_range = pd.date_range(ffwsjlow.index.min(), sofr.index.max(), freq="D")
# each series on the daily index, forward-filled (weekends, holidays, monthly values)
s_wsj = ffwsjlow["ffwsjlow"].reindex(full_date_range).ffill()
s_fed = fed_funds["FED eff. rate"].reindex(full_date_range).ffill() # NaN before its first observation
s_sofr = sofr["SOFR"].reindex(full_date_range).ffill() # NaN before its first observation
# pick the source by period: WSJ -> Fed funds -> SOFR
ff_start, sofr_start = fed_funds.index.min(), sofr.index.min()
bc_yearly = s_wsj.where(full_date_range < ff_start,
s_fed.where(full_date_range < sofr_start, s_sofr))
borrowing_cost = pd.DataFrame({
"ffwsjlow": s_wsj.where(full_date_range < ff_start),
"FED eff. rate": s_fed.where(full_date_range < sofr_start),
"SOFR": s_sofr,
"bc_yearly": bc_yearly,
})
borrowing_cost.to_csv(r"output\borrowing_cost.csv")
# daily rate for the leverage model (rates are in percent)
borrowing_cost["bc"] = (1 + borrowing_cost["bc_yearly"] / 100) ** (1 / 365) - 1
borrowing_cost[["ffwsjlow", "FED eff. rate", "SOFR"]].plot(title="Federal Funds Rate and SOFR", ylabel="Rate (%)", xlabel="Date", figsize=(8, 5))
plt.figure()
borrowing_cost["bc_yearly"].plot(title = "Borrowing Cost in basis points p.a.", figsize=(8,5))
data_full = returns[["return_daily"]].merge(borrowing_cost["bc"], left_index = True, right_index = True, how = "outer")
data_full["er"] = 1.0089**(1/365) - 1
# data_full["a"] = 0 # calculate later
data_full["a"] = (1.011911)**(1/365) - 1
data_full.to_csv(r"output\SP500_basis_leveraged_returns.csv")
In [753]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter, PercentFormatter
## Functions to Plot leveraged returns for arbitrary start value and start/end-date
def cummulate_leveraged_returns(
returns_df,
start_value=10000.0,
leverage = 3,
start_date = pd.to_datetime("1932-06-01"), end_date= pd.to_datetime("2026-10-01"),
returns = "return_daily", expense_ratio = "er", borrowing_cost = "bc", adjustment_factor = "a",
output = False, output_name = "results.xlsx"
):
def calc_leveraged_returns(L, r, er, bc, a):
return L*r - er - (L-1)*bc - a
dates_mask = (
returns_df[returns].notna()
& returns_df[expense_ratio].notna()
& returns_df[borrowing_cost].notna()
& returns_df[adjustment_factor].notna()
& (returns_df.index >= start_date)
& (returns_df.index <= end_date)
)
returns_available = returns_df[dates_mask]
returns_available.loc[returns_df[dates_mask].index.min(), returns] = 0
cum_og_returns = (1 + returns_available.loc[dates_mask, returns]).cumprod()
max_dd_og = (cum_og_returns / cum_og_returns.cummax() - 1).min()
leveraged_returns = calc_leveraged_returns(L = leverage, r = returns_available[returns], er = returns_available[expense_ratio], bc = returns_available[borrowing_cost], a = returns_available[adjustment_factor])
# set starting value to 0 such tha cummulated returns start at 1
leveraged_returns.iloc[0] = 0
cum_leveraged_returns = (1 + leveraged_returns).cumprod()
er_annual = (1 + returns_available[expense_ratio].dropna().iloc[0])**365 - 1
adj_factor = (1 + returns_available[adjustment_factor].dropna().iloc[0])**365 - 1
# max percentage loss from peak
max_dd = (cum_leveraged_returns / cum_leveraged_returns.cummax() - 1).min()
if output:
# export results
with pd.ExcelWriter(rf"output/{output_name}", engine="openpyxl") as writer:
res = pd.DataFrame({
"lev_returns": leveraged_returns,
"cum_lev_returns": cum_leveraged_returns * start_value,
"er_annual": returns_available[expense_ratio],
"adjustment_factor": returns_available[adjustment_factor],
"bc": returns_available[borrowing_cost]
})
res.to_excel(writer, sheet_name=f"data", index_label="date")
pd.DataFrame({
"leverage": [leverage],
"max_drawdown": [max_dd],
"og_max_drawdown": [max_dd_og],
"start_value": [start_value],
"start_date": [start_date],
"end_date": [end_date],
}, index=["value"]).to_excel(writer, sheet_name="metadata")
return {
"lev_returns": leveraged_returns,
"cum_lev_returns": cum_leveraged_returns * start_value,
"er_annual": returns_available[expense_ratio],
"adjustment_factor": adj_factor,
"max_drawdown": max_dd,
"start_value": start_value,
"cum_og_returns": cum_og_returns * start_value,
"og_max_drawdown": max_dd_og,
"bc": returns_available[borrowing_cost],
}
def human(v, _):
for div, suffix in ((1e9, "B"), (1e6, "M"), (1e3, "k")):
if abs(v) >= div:
return f"{v / div:g}{suffix}"
return f"{v:g}"
def plot_leveraged_returns(
cum_lev_returns, er_annual, adjustment_factor, max_drawdown=None, start_value=1.0,
title="SP500 leveraged returns", x_label="Date", y_label="Accumulated Returns",
):
cum = cum_lev_returns.to_frame("cum") if isinstance(cum_lev_returns, pd.Series) else cum_lev_returns
show_dd = max_drawdown is not None
if show_dd:
fig, (ax, ax_dd) = plt.subplots(
2, 1, figsize=(12, 8), sharex=True,
gridspec_kw={"height_ratios": [2, 1], "hspace": 0.05},
)
else:
fig, ax = plt.subplots(figsize=(12, 6))
ax_dd = None
# top: wealth on a log axis
for col in cum.columns:
ax.plot(cum.index, cum[col], label=str(col))
colors = [line.get_color() for line in ax.get_lines()]
ax.set_yscale("log")
ax.yaxis.set_major_formatter(FuncFormatter(human))
ax.axhline(start_value, color="gray", lw=0.8, ls="--")
ax.set_ylabel(y_label)
ax.set_title(title)
ax.grid(True, which="major", alpha=0.3)
if cum.shape[1] > 1:
ax.legend()
footer = f"ER = {er_annual:.2%} Adjustment Factor = {adjustment_factor:.2%}"
if show_dd:
dd = cum / cum.cummax() - 1 # underwater curve, <= 0
for col, color in zip(dd.columns, colors):
ax_dd.plot(dd.index, dd[col], color=color, lw=0.8)
ax_dd.fill_between(dd.index, dd[col], 0, color=color, alpha=0.15)
ax_dd.set_ylim(dd.min().min() * 1.05, 0)
ax_dd.yaxis.set_major_formatter(PercentFormatter(1.0, decimals=0))
ax_dd.set_ylabel("Drawdown")
ax_dd.set_xlabel(x_label)
ax_dd.grid(True, alpha=0.3)
# scalar -> "max drawdown = x"; several values -> all of them, labeled if keyed
if isinstance(max_drawdown, (dict, pd.Series)):
items = dict(max_drawdown).items()
dd_text = ", ".join(f"{str(k)}: {-abs(v):.2%}" for k, v in items)
elif np.ndim(max_drawdown) > 0:
dd_text = ", ".join(f"{-abs(v):.2%}" for v in max_drawdown)
else:
dd_text = f"{-abs(max_drawdown):.2%}"
footer += f" max drawdown = {dd_text}"
else:
ax.set_xlabel(x_label)
fig.subplots_adjust(bottom=0.12)
fig.text(0.5, 0.02, footer, ha="center")
return ax, ax_dd
2.3 Visualizing and drawdowns¶
In order to visualize the extreme drawdowns in the leverage, we use the maximum drawdown. Let $S_t$ be the value of cummulated returns at day $t$ and $M_t$ the maximum up to time $t$. Let $D_t$ be the percentage loss from the peak.
$$ M_t = \max_{u\leq t}S_u \ \ \ \ D_t = \frac{S_t}{M_t}-1$$
We define the max drawdown as most negative percentage loss from last peak at time t.
$$ \text{MaxDD} = \min_t{D_t}$$
In [765]:
data_full = pd.read_csv(r"output\SP500_basis_leveraged_returns.csv", index_col = 0)
data_full.index = pd.to_datetime(data_full.index)
# calc leverage
leverages = [1.5, 2, 3]
results = {L: cummulate_leveraged_returns(data_full, leverage=L) for L in leverages}
cum_all = pd.DataFrame({f"{L}x": r["cum_lev_returns"] for L, r in results.items()})
ret_all = pd.DataFrame({f"{L}x": r["lev_returns"] for L, r in results.items()})
cum_all.insert(0, "og 1x", results.get(leverages[0]).get("cum_og_returns"))
first = results[leverages[0]] # er, adjustment factor and start value are the same for all leverages
drawdowns = {
"og 1x": first.get("og_max_drawdown"),
**{key: value["max_drawdown"] for key, value in results.items()},
}
# export results
# with pd.ExcelWriter(r"output/SP500_leverage_results.xlsx", engine="openpyxl") as writer:
# data_full.to_excel(writer, sheet_name="basis data", index_label="date")
# for L, r in results.items():
# res = pd.DataFrame({
# "lev_returns": r["lev_returns"],
# "cum_lev_returns": r["cum_lev_returns"],
# })
# res.to_excel(writer, sheet_name=f"SP500 x{L}", index_label="date")
# plot results
ax, ax_dd = plot_leveraged_returns(
cum_lev_returns=cum_all,
er_annual=first["er_annual"].iloc[0],
adjustment_factor=first["adjustment_factor"],
max_drawdown=drawdowns,
start_value=first["start_value"],
)
2.4 Validate 3x Leveraged Returns against UPRO¶
In order to check if calculated returns are suitable, we check the 3x SP500 against UPRO. We have UPRO data since 2008. We set the adjustment factor to be the avarage difference between synthetical 3x returns (with adjustment factor 0) and UPRO returns.
2.4.1 Approximate adjustment factor¶
The constant annual adjustment $a^*$ is found numerically by minimising the sum of squared differences between the synthetic 3x wealth path $S_t(a)$ and the UPRO path $S^{\text{UPRO}}_t$, i.e.
$$a^* = \arg\min_a \sum_t \big(S_t(a) - S^{\text{UPRO}}_t\big)^2$$
In [766]:
from scipy.optimize import minimize_scalar
start_value = 10000
L = 3
# get UPRO data
up = upro.set_index("Date").rename(columns={"Close":"UPRO"})
up["upro_returns"] = up["UPRO"].pct_change()
up["UPRO_st"] = (up["UPRO"] / up["UPRO"].iloc[0]) * start_value
# get SP500
data_full = pd.read_csv(r"output\SP500_basis_leveraged_returns.csv", index_col = 0)
data_full.index = pd.to_datetime(data_full.index)
# common window and the a = 0 synthetic return
d = data_full.join(up[["UPRO_st"]], how="inner")
d = d.dropna(subset=["return_daily", "bc", "er", "UPRO_st"])
r0 = (L * d["return_daily"] - d["er"] - (L - 1) * d["bc"]).to_numpy()
# UPRO price
upro_st = d["UPRO_st"].to_numpy()
# calculate adjustment factor to minimize sse between upro price and cumulative returns
def sse(a_annual):
a_row = (1 + a_annual) ** (1 / 365) - 1
synth = start_value * np.cumprod(1 + r0 - a_row)
return np.sum((synth - upro_st) ** 2) / start_value**2 # scaled to avoid huge numbers
res = minimize_scalar(sse, bounds=(-0.05, 0.20), method="bounded", options={"xatol": 1e-10})
a_opt = res.x
print(f"optimal annual adjustment: {a_opt:.4%} (SSE = {res.fun:.4g})")
optimal annual adjustment: 1.1911% (SSE = 1.406e+04)
2.4.2 Compare synthetic 3x vs UPRO¶
In [768]:
# get UPRO data
up = upro.set_index("Date").rename(columns={"Close":"UPRO"})
up["upro_returns"] = up["UPRO"].pct_change()
up["UPRO_st"] = (up["UPRO"] / up["UPRO"].iloc[0]) * start_value
# get SP500
data_full = pd.read_csv(r"output\SP500_basis_leveraged_returns.csv", index_col = 0)
data_full.index = pd.to_datetime(data_full.index)
d = data_full.join(up[["UPRO_st", "upro_returns"]])
# plot UPRO against SP500 with calculated adjustment factor
lev_returns_dict = cummulate_leveraged_returns(d, leverage = 3, start_value = start_value, start_date = up.index[0], end_date = up.index[-1])
cum_lev_returns = lev_returns_dict.get("cum_lev_returns")
# Plot snyth vs UPRO cum-returns
cum_all = pd.concat([cum_lev_returns, up["UPRO_st"]], axis = 1).rename(columns={0: "synth 3x", "UPRO_st": "UPRO"})
ax, ad = plot_leveraged_returns(
cum_lev_returns=cum_all,
er_annual=lev_returns_dict["er_annual"].iloc[0],
adjustment_factor=lev_returns_dict["adjustment_factor"],
max_drawdown=None, # worst case, shown in the footer
start_value=lev_returns_dict["start_value"],
title = "SP500 3x vs UPRO"
)
# Plot return differences for
delta_returns = d["return_daily"] - d["upro_returns"]
q99 = delta_returns.quantile(0.99)
q01 = delta_returns.quantile(0.01)
filtered = delta_returns[(delta_returns <= q99) & (delta_returns > q01)]
fig, (ax_density, ax_returns) = plt.subplots(1, 2, figsize=(14, 5))
filtered.plot(
kind="density",
ax=ax_density,
title="Density of Daily Return Differences (UPRO vs. S&P 500 3x)\n1st-99th Percentile",
xlabel="Daily return difference (decimal)",
ylabel="Density",
)
filtered.plot(
ax=ax_returns,
title="Daily Return Differences\n1st-99th Percentile",
xlabel="Date",
ylabel="Difference",
)
filtered_mean = filtered.mean()
filtered_std = filtered.std()
ax_returns.axhline(filtered_mean,color="red",linestyle="--",label=f"Mean = {filtered_mean:.6f}, Std = {filtered_std:.6f}",)
ax_returns.legend()
fig.tight_layout()
plt.show()
3. Result Analysis¶
3.1. Single investment for different holding periods¶
Assume that for every year in the period (1952-2026) two equally big investments in the 3x SP500 and non-leveraged SP500 are taken on a random day $t=0$ in that year. Assume further, that this investment is hold for $H\in \{1,10,20,30\}$ years. Compare the returns of the investments by dividing leveraged-returns by non-leveraged returns. Let $P(t)$ be the cummulated percentage return of the investments. Also compar the max drawdown of investments in that timeframe.
$$ r_{\text{total}} = \frac{P(H)}{P(0)}, \ \ \ \ r_{\text{ratio}}=\frac{r_{\text{total 3x}}}{r_{\text{total 1x}}}$$
The followig
In [ ]:
data_full = pd.read_csv(r"output\SP500_basis_leveraged_returns.csv", index_col = 0)
data_full.index = pd.to_datetime(data_full.index)
d = data_full.dropna(subset=["return_daily", "bc", "er"])
H_list = [1, 10, 20, 30]
rng = np.random.default_rng(42) # change the seed for a different sequence
first_date, last_date = d.index.min(), d.index.max()
years = d.index.year.unique() # all calendar years in the data
results = {}
for H in H_list:
# set start_dates randomly: one investment_date per year
latest = last_date - pd.DateOffset(years=H) # latest start date so that the window still fits
start_dates = []
for year in years:
lo = max(pd.Timestamp(year=year, month=1, day=1), first_date)
hi = min(pd.Timestamp(year=year, month=12, day=31), latest)
if lo > hi: # no valid start date in this year for this H
continue
offset = int(rng.integers((hi - lo).days + 1)) # random day between lo and hi (inclusive)
start_dates.append((lo + pd.Timedelta(days=offset)).normalize())
start_dates = pd.DatetimeIndex(start_dates) # investment_date for each year
# get results
rows = []
for s in start_dates:
r = cummulate_leveraged_returns(d, leverage=L, start_date=s, end_date=s + pd.DateOffset(years=H))
rows.append({
"start": s,
"end": s + pd.DateOffset(years=H),
"lev_mult": r["cum_lev_returns"].iloc[-1] / r["start_value"],
"og_mult": r["cum_og_returns"].iloc[-1] / r["start_value"],
"lev_dd": r["max_drawdown"],
"og_dd": r["og_max_drawdown"],
})
res = pd.DataFrame(rows).set_index("start")
res["ratio"] = res["lev_mult"] / res["og_mult"]
results[H] = res
# import openpyxl
# with pd.ExcelWriter(r"output\holding_period_results.xlsx", engine="openpyxl") as writer:
# for H, res in results.items():
# res.to_excel(writer, sheet_name=f"{H} year", index_label="start")
In [770]:
for k in range(0,len(H_list)):
fig, axs = plt.subplots(1, 2, figsize=(13, 5))
fig.suptitle(f"{L}x Leverage vs. 1x: Return Ratio and Drawdown for Holding Period of {H_list[k]}Y")
res = results.get(H_list[k])
# b) leverage vs unleveraged: where does leverage win?
ax = axs[0]
ax.plot(res.index, res["ratio"], color="black", lw=0.8)
ax.axhline(1, color="gray", ls="--")
ax.fill_between(res.index, res["ratio"], 1, where=res["ratio"] < 1, color="red", alpha=0.3, label="leverage loses")
ax.fill_between(res.index, res["ratio"], 1, where=res["ratio"] >= 1, color="green", alpha=0.2, label="leverage wins")
avg_ratio = res["ratio"].mean()
ax.axhline(avg_ratio, color="blue", ls=":", label=f"mean ratio = {avg_ratio:.2f}")
ax.set_yscale("log")
ax.set_title(f"(Total Return {L}x) / (Total Return 1x)"); ax.legend()
ax.set_xlabel("year of investment")
ax.set_ylabel("log ratio (3x / 1x)")
# c) worst drawdown inside each window
ax = axs[1]
ax.plot(res.index, res["lev_dd"], label=f"{L}x")
ax.plot(res.index, res["og_dd"], label="1x")
ax.yaxis.set_major_formatter(plt.matplotlib.ticker.PercentFormatter(1.0, decimals=0))
ax.set_title(f"Max Drawdown within Investment-Window"); ax.legend()
ax.set_xlabel("year of investment")
plt.tight_layout()
plt.show()
3.2 Invest every month¶
Assume that 200€ are invested at at a random day in a month for 10 years or 20 years in the 3xSP500 and SP500. As a starting period, choose the following dates:
| Start | Window | Why it was tough | | 1966-01 | 1966–1976 | Start near the 1966 peak. Nominal index went sideways for years, with several bear markets (1966, 1968–70, 1973–74). Inflation and interest rates kept rising, so borrowing costs ate the leverage. | | 1969-01 | 1969–1979 | Start just after the 1968 peak. Bear market 1969–70, then the 1973–74 crash (about -48%), followed by stagflation. Weak real returns and high rates. | | 1973-01 | 1973–1983 | Start right at the top before the 1973–74 crash (about -48%), then oil shocks and Volcker-era rates (Fed funds near 19% in 1981). The worst case for a leveraged strategy. | | 2000-09 | 2000–2010 | Start near the dotcom peak. Two bear markets (2000–02 and 2007–09, each around -50%) and a 10-year return close to zero ("lost decade"). |
In [771]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter
L, H, MONTHLY = 3, 10, 200.0
MIN_GAP = 6
def dca_curves(df, start, H = H, L = L, monthly = MONTHLY, seed = 0):
"""Sparplan mit zufälligem Einzahlungstag pro Monat: Wert (3x, 1x) und Gesamteinzahlung."""
end = start + pd.DateOffset(years=H)
res = cummulate_leveraged_returns(df, start_value=1.0, leverage=L, start_date=start, end_date=end)
G_lev, G_og = res["cum_lev_returns"], res["cum_og_returns"] # beide starten bei 1
idx = G_lev.index
rng = np.random.default_rng(seed) # gleicher Seed -> gleiche Einzahlungstage
ms = pd.date_range(start, periods=H * 12 + 1, freq="MS")
offs = (rng.random(H * 12) * np.asarray(ms[:-1].days_in_month)).astype(int)
c_dates = ms[:-1] + pd.to_timedelta(offs, unit="D") # zufälliger Tag in jedem Monat
pos = np.minimum(idx.searchsorted(c_dates), len(idx) - 1) # erster verfügbarer Handelstag >= Einzahlungstag
def value(G):
g = G.to_numpy()
shares = np.zeros(len(idx))
np.add.at(shares, pos, monthly / g[pos]) # gekaufte Anteile je Einzahlung
return pd.Series(g * shares.cumsum(), index=idx)
invested = pd.Series(monthly * np.cumsum(np.bincount(pos, minlength=len(idx))), index=idx)
return value(G_lev), value(G_og), invested
starts = ["1966-01-01", "1969-01-01", "1973-01-01", "2000-09-01"]
euro = FuncFormatter(lambda v, _: f"{v:,.0f} €")
fig, axs = plt.subplots(2, 2, figsize=(13, 7))
summary = []
for ax, s in zip(axs.ravel(), starts):
s = pd.Timestamp(s)
v_lev, v_og, inv = dca_curves(data_full, s, H=H, L=L, monthly=MONTHLY, seed=s.year)
ax.plot(v_lev, label=f"{L}x S&P 500: {v_lev.iloc[-1]:,.0f} €")
ax.plot(v_og, label=f"S&P 500: {v_og.iloc[-1]:,.0f} €")
ax.plot(inv, "k--", label=f"Eingezahlt: {inv.iloc[-1]:,.0f} €")
ax.yaxis.set_major_formatter(euro)
ax.set_title(f"{s:%Y-%m} bis {s + pd.DateOffset(years=H):%Y-%m}")
ax.grid(True, alpha=0.3)
ax.legend(loc="upper left")
summary.append({
"Start": f"{s:%Y-%m}",
"3x": v_lev.iloc[-1], "1x": v_og.iloc[-1], "eingezahlt": inv.iloc[-1],
"3x/eingezahlt": v_lev.iloc[-1] / inv.iloc[-1],
"1x/eingezahlt": v_og.iloc[-1] / inv.iloc[-1],
"3x min. Wert/eingezahlt": (v_lev / inv).min(), # tiefster Punkt relativ zum bis dahin Eingezahlten
})
fig.suptitle(f"{MONTHLY:.0f} € pro Monat (zufälliger Tag), {H} Jahre: {L}x S&P 500 vs. S&P 500", y=1.0)
plt.tight_layout(); plt.show()
In [772]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter
L, H, MONTHLY = 3, 20, 200.0
MIN_GAP = 6
starts = ["1966-01-01", "1969-01-01", "1973-01-01", "2000-09-01"]
euro = FuncFormatter(lambda v, _: f"{v:,.0f} €")
fig, axs = plt.subplots(2, 2, figsize=(13, 7))
summary = []
for ax, s in zip(axs.ravel(), starts):
s = pd.Timestamp(s)
v_lev, v_og, inv = dca_curves(data_full, s, H=H, L=L, monthly=MONTHLY, seed=s.year)
ax.plot(v_lev, label=f"{L}x S&P 500: {v_lev.iloc[-1]:,.0f} €")
ax.plot(v_og, label=f"S&P 500: {v_og.iloc[-1]:,.0f} €")
ax.plot(inv, "k--", label=f"Eingezahlt: {inv.iloc[-1]:,.0f} €")
ax.yaxis.set_major_formatter(euro)
ax.set_title(f"{s:%Y-%m} bis {s + pd.DateOffset(years=H):%Y-%m}")
ax.grid(True, alpha=0.3)
ax.legend(loc="upper left")
summary.append({
"Start": f"{s:%Y-%m}",
"3x": v_lev.iloc[-1], "1x": v_og.iloc[-1], "eingezahlt": inv.iloc[-1],
"3x/eingezahlt": v_lev.iloc[-1] / inv.iloc[-1],
"1x/eingezahlt": v_og.iloc[-1] / inv.iloc[-1],
"3x min. Wert/eingezahlt": (v_lev / inv).min(), # tiefster Punkt relativ zum bis dahin Eingezahlten
})
fig.suptitle(f"{MONTHLY:.0f} € pro Monat (zufälliger Tag), {H} Jahre: {L}x S&P 500 vs. S&P 500", y=1.0)
plt.tight_layout(); plt.show()
4. 3x 30Y Treasury ETF¶
A young retail investor suggested in 2019 to hedge the huge downside risk of 3x SP500 with 3x leverage on 30Y Treasury Yield fund TYX (the post is linked in [6]). This fund was set up in 1985 [7]. Data for Market Yield on U.S. Treasury Securities at 30-Year Constant Maturity at FED starts in 1977 ([8]). A user in the blog noted above fortunately typed monthly treasury yields starting from 1941 given in a book into a spreadsheet ([9]).
We interpolate data from 1941 until 1977 from monthly to daily by linear interpolation.
$$ y(t) = y(t_0) + \frac{t-t_0}{t_1-t_0}\left(y(t_1)-y(t_0)\right), \ \ \ \text{for all }t\in (t_0,t_1) $$
In [773]:
yields = pd.read_csv(r"data\yields_eom_1941.csv", sep = ";", decimal = ",")
yields["End of Month"] = (pd.to_datetime(yields["End of Month"].str.replace(r"-(\d{2})$", r"-19\1", regex=True),format="%b-%Y",) + pd.offsets.MonthEnd(0))
yields = yields.set_index(yields["End of Month"]).drop(columns=["End of Month", "Unnamed: 12"]).dropna()
# interpolate to daily returns
full_date_range = pd.date_range(yields.index.min(), yields.index.max(), freq="D")
yields = yields.reindex(full_date_range).interpolate(method="time")
y30 = yields["Long"]
# FED data is given daily. Fill missing values by copying the last value
dgs30 = read_fed(r"data\DGS30.csv", date_col = "observation_date", value_col = "value").set_index("observation_date")
full_date_range = pd.date_range(dgs30.index.min(), dgs30.index.max(), freq="D")
dgs30 = dgs30.reindex(full_date_range).ffill()
# combine datasets
dat = dgs30.join(yields[["Long"]], how = "outer").rename(columns = {
"Long": "30Y Monthly Treasury Yields (book)",
"value": "30Y Daily Treasury Yields (FED)"
})
date_range_1 = pd.date_range(yields.index.min(), dgs30.index.min(), freq="D")
date_range_2 = pd.date_range(dgs30.index.min(), dgs30.index.max(), freq="D")
dat.loc[date_range_1, "30Y"] = dat.loc[date_range_1, "30Y Monthly Treasury Yields (book)"]
dat.loc[date_range_2, "30Y"] = dat.loc[date_range_2, "30Y Daily Treasury Yields (FED)"]
dat.to_csv(r"output\historical_30Y_treasury_yields.csv")
# plot results
fig, axs = plt.subplots(1, 2, figsize=(13, 5))
axs[0].plot(dat[["30Y Daily Treasury Yields (FED)", "30Y Monthly Treasury Yields (book)"]])
axs[0].set_title("30Y Interpolated Book Yields vs Daily FED yields")
diff = (dat["30Y Monthly Treasury Yields (book)"] - dat["30Y Daily Treasury Yields (FED)"]).dropna()
diff_qt = diff[(diff <= diff.quantile(0.99)) & (diff >= diff.quantile(0.01))]
diff_mean = diff.mean()
axs[1].plot(diff_qt)
axs[1].axhline(diff_mean, color="blue", ls=":", label=f"mean = {diff_mean:.4f}, std deviation = {diff.std():.4f}")
axs[1].legend()
axs[1].set_title("monthly interpolated to daily (book) minus daily (FED)\n(1st–99th Percentiles)")
Out[773]:
Text(0.5, 1.0, 'monthly interpolated to daily (book) minus daily (FED)\n(1st–99th Percentiles)')
In [774]:
start_value = 10000
# Goal: Hit the 20Y+ Bond ETF (TYX)
tyx = read_yfinance(r"data/TYX.csv").set_index("Date").rename(columns={"Close":"TYX"})
full_date_range = pd.date_range(tyx.index.min(), tyx.index.max(), freq="D")
tyx = tyx.reindex(full_date_range).ffill()
tyx["tyx_returns"] = tyx["TYX"].pct_change()
tyx.loc[tyx.index.min(), "tyx_returns"] = 0
tyx["tyx_cum_returns"] = (1+tyx["tyx_returns"]).cumprod() * start_value
#
hist_30y = pd.read_csv(r"output\historical_30Y_treasury_yields.csv", index_col = 0)
hist_30y.index = pd.to_datetime(hist_30y.index)
hist_30y.loc[tyx.index.min():,["30Y"]] = (hist_30y.loc[tyx.index.min():,["30Y"]] / hist_30y.loc[tyx.index.min(),["30Y"]]) * start_value
hist_30y["30Y_returns"] = hist_30y["30Y"].pct_change()
# plot results
fig, axs = plt.subplots(1, 1, figsize=(6, 5))
dat = tyx[["tyx_cum_returns"]].join(hist_30y[["30Y"]], how = "inner")
axs.plot(dat)
axs.set_title("30Y Interpolated Book Yields vs 20Y+ Yield Index")
axs.legend()
# dat1 = hist_30y[["30Y_returns"]].join(tyx["tyx_returns"], how = "inner")
# diff = (dat1["30Y_returns"] - dat1["tyx_returns"])
# # diff_mean = diff.mean()
# axs[1].plot(diff_qt)
# # axs[1].axhline(diff_mean, color="blue", ls=":", label=f"mean = {diff_mean:.4f}, std deviation = {diff.std():.4f}")
# # axs[1].legend()
# axs[1].set_title("monthly interpolated to daily (book) minus daily index returns (TYX)\n(1st–99th Percentiles)")
# diff
C:\Users\fabia\AppData\Local\Temp\ipykernel_11580\2069499147.py:22: UserWarning: No artists with labels found to put in legend. Note that artists whose label start with an underscore are ignored when legend() is called with no argument. axs.legend()
Out[774]:
<matplotlib.legend.Legend at 0x21299aacc20>
To be continued¶
Next goal is to derive returns of treasury funds using yield curve