Implement the Multivariate Bootstrap method used in Aizhan Anarkulova's articles [INCOMPLETE]¶
This is "101" implementation, playing with simpler data sets and fewer assets to understand the power and the limitations of the bootstrap method.
- Respects serial dependence and cross-asset correlation via synchronized stationary block bootstrap
- Downloads free adjusted-close price data via yfinance
- Generates grid of allocations (10% increments) over SPY, VT, GLD, IEF
- Estimates expected outcomes over 5, 10, 20, and 30-year horizons: mean final value, median final value, 5th/95th percentiles, probability of loss, mean/median CAGR
- Automatically handles simulation start point based on overlapping data availability
In [1]:
import numpy as np
import pandas as pd
from dataclasses import dataclass
from typing import List, Dict
from tqdm import tqdm
import plotly.express as px
import yfinance as yf
yf.__version__
Out[1]:
In [7]:
@dataclass
class BootstrapOutcome:
horizon_years: int
weights: Dict[str, float]
weights_txt: str # For display purposes
mean_final: float
median_final: float
p5_final: float
p95_final: float
probability_loss: float
mean_cagr: float
median_cagr: float
def download_adjusted_prices(tickers: List[str], start: str = None, end: str = None, interval: str = "1d") -> pd.DataFrame:
data = yf.download(tickers, start=start, end=end, interval=interval, progress=False, auto_adjust=True)
if 'Close' in data.columns:
prices = data['Close'].copy()
else:
prices = data.copy()
if isinstance(prices, pd.Series):
prices = prices.to_frame()
prices.columns = tickers if len(tickers) > 1 else [tickers[0]]
prices = prices.sort_index()
return prices
def compute_monthly_returns(prices: pd.DataFrame) -> pd.DataFrame:
monthly = prices.resample('ME').last() # ME = Month End
returns = monthly.pct_change().dropna(how='any')
return returns
def stationary_bootstrap_indices(n_obs: int, horizon: int, avg_block_size: int, rng: np.random.Generator) -> np.ndarray:
p_new = 1.0 / avg_block_size
indices = np.empty(horizon, dtype=int)
indices[0] = rng.integers(0, n_obs)
for t in range(1, horizon):
if rng.random() < (1 - p_new):
indices[t] = (indices[t-1] + 1) % n_obs
else:
indices[t] = rng.integers(0, n_obs)
return indices
def simulate_bootstrap_paths(returns: pd.DataFrame, horizon_months: int, n_sims: int,
avg_block_size: int, rng: np.random.Generator) -> np.ndarray:
data = returns.to_numpy()
n_obs, n_assets = data.shape
block = min(max(2, avg_block_size), max(2, n_obs))
sims = np.empty((n_sims, horizon_months, n_assets))
for i in range(n_sims):
idx = stationary_bootstrap_indices(n_obs, horizon_months, block, rng)
sims[i] = data[idx]
return sims
def generate_weight_grid(n_assets: int, step: float = 0.10) -> List[np.ndarray]:
from itertools import product
steps = np.arange(0, 1 + step, step)
combinations = product(steps, repeat=n_assets)
grids = []
for combo in combinations:
if abs(sum(combo) - 1.0) < 1e-10:
grids.append(np.array(combo))
return grids
def compute_portfolio_statistics(sims: np.ndarray, weights: np.ndarray, years: int) -> BootstrapOutcome:
"""
sims: shape (n_sims, n_months, n_assets)
weights: shape (n_assets,)
years: investment horizon in years (n_months = years * 12)
Returns all statistics based on annualized (CAGR) returns.
"""
sims_pr = sims @ weights # shape (n_sims, n_months)
# Compute final value for each simulation
final_values = np.prod(1 + sims_pr, axis=1)
# Compute CAGR for each simulation
cagr = final_values ** (1.0 / years) - 1.0
# All statistics below are based on annualized returns (CAGR)
return BootstrapOutcome(
horizon_years=years,
weights=dict(zip(['SPY', 'VT', 'GLD', 'IEF'], weights)),
weights_txt=' + '.join(f"{w:.0%} {k}" for w, k in zip(weights, ['SPY', 'VT', 'GLD', 'IEF']) if w > 0),
mean_final=float(np.mean(cagr)),
median_final=float(np.median(cagr)),
p5_final=float(np.percentile(cagr, 5)),
p95_final=float(np.percentile(cagr, 95)),
probability_loss=float(np.mean(cagr < 0)),
mean_cagr=float(np.mean(cagr)),
median_cagr=float(np.median(cagr)),
)
In [10]:
tickers = ['SPY', 'VT', 'GLD', 'IEF']
start = start="1900-01-01"
end = None
n_sims = 4000
avg_block_months = 12
horizons = [5, 10, 20, 30]
weight_step = 0.10
random_seed = 42
rng = np.random.default_rng(random_seed)
prices = download_adjusted_prices(tickers, start, end)
print(f"Display prices: oldest available since {start}")
display(prices)
returns = compute_monthly_returns(prices)
print(f"Display returns: oldest available FOR ALL since {start}")
display(returns)
if returns.shape[0] < 36:
raise ValueError(f"Not enough data points: only {returns.shape[0]} months available.")
effective_block = min(max(2, avg_block_months), max(2, returns.shape[0] // 5))
print(f"Monthly returns from {returns.index.min().date()} to {returns.index.max().date()} — {returns.shape[0]} months")
print(f"Effective average block size (months): {effective_block}")
weight_sets = generate_weight_grid(n_assets=len(tickers), step=weight_step)
print(f"Generated {len(weight_sets)} weight sets with step {weight_step:.2f} for {len(tickers)} assets")
outcomes = []
total = len(horizons) * len(weight_sets)
with tqdm(total=total, desc="Simulations") as pbar:
for years in horizons:
months = years * 12
sims = simulate_bootstrap_paths(returns, months, n_sims, effective_block, rng)
for w in weight_sets:
outcome = compute_portfolio_statistics(sims, w, years)
outcomes.append(outcome)
pbar.update(1)
df = pd.DataFrame([{
'horizon_years': outcome.horizon_years,
'weight_txt': outcome.weights_txt,
'mean_final': outcome.mean_final,
'median_final': outcome.median_final,
'p5_final': outcome.p5_final,
'p95_final': outcome.p95_final,
'probability_loss': outcome.probability_loss,
'mean_cagr': outcome.mean_cagr,
'median_cagr': outcome.median_cagr,
**{f'w_{ticker}': outcome.weights[ticker] for ticker in tickers}
} for outcome in outcomes])
df = df.round(3).sort_values(['horizon_years', 'mean_cagr'], ascending=[True, False]).reset_index(drop=True)
display(df)
In [11]:
for ww in ['w_SPY', 'w_VT', 'w_GLD', 'w_IEF']:
fig = px.scatter(
df,
x="median_final",
y="p5_final",
color=ww,
facet_col="horizon_years",
hover_data=[col for col in df.columns if not col.startswith('w_') and col not in ['median_final', 'p5_final', 'horizon_years']],
)
fig.update_layout(height=250, width=1200, margin=dict(t=20, b=20, l=20, r=20))
fig.for_each_yaxis(lambda y: y.update(matches=None, showticklabels=True, tickformat=".0%"))
fig.for_each_xaxis(lambda x: x.update(matches=None, tickformat=".0%"))
fig.show()