import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
class MonteCarloSimulation:
def __init__(self,
initial_balance: float,
annual_contribution: float,
contribution_growth_rate: float,
retirement_year: int,
annual_withdrawal: float,
withdrawal_growth_rate: float,
years: int,
num_simulations: int,
inflation_rate: float,
tax_rate: float,
asset_allocation: dict,
use_bootstrap: bool = False,
historical_returns: pd.DataFrame = None):
"""
asset_allocation: {'stocks': 0.7, 'bonds': 0.3}
historical_returns: DataFrame with 'stocks' and 'bonds' columns if bootstrapping
"""
self.initial_balance = initial_balance
self.annual_contribution = annual_contribution
self.contribution_growth_rate = contribution_growth_rate
self.retirement_year = retirement_year
self.annual_withdrawal = annual_withdrawal
self.withdrawal_growth_rate = withdrawal_growth_rate
self.years = years
self.num_simulations = num_simulations
self.inflation_rate = inflation_rate
self.tax_rate = tax_rate
self.asset_allocation = asset_allocation
self.use_bootstrap = use_bootstrap
self.historical_returns = historical_returns
# Defaults if bootstrapping not used
self.asset_return_means = {'stocks': 0.07, 'bonds': 0.03}
self.asset_return_stds = {'stocks': 0.15, 'bonds': 0.05}
def simulate_return(self):
if self.use_bootstrap and self.historical_returns is not None:
sample = self.historical_returns.sample(n=1, replace=True).iloc[0]
portfolio_return = sum(sample[asset] * weight for asset, weight in self.asset_allocation.items())
else:
portfolio_return = 0.0
for asset, weight in self.asset_allocation.items():
mean = self.asset_return_means[asset]
std = self.asset_return_stds[asset]
asset_return = np.random.normal(mean, std)
portfolio_return += asset_return * weight
return portfolio_return
def run(self):
simulations = []
for sim in range(self.num_simulations):
balance = self.initial_balance
yearly_balances = []
for year in range(self.years):
annual_return = self.simulate_return()
balance *= (1 + annual_return)
if year < self.retirement_year:
contribution = self.annual_contribution * ((1 + self.contribution_growth_rate) ** year)
balance += contribution
else:
withdrawal = self.annual_withdrawal * ((1 + self.withdrawal_growth_rate) ** (year - self.retirement_year))
taxed_withdrawal = withdrawal / (1 - self.tax_rate)
balance -= taxed_withdrawal
# Inflation adjustment
balance /= (1 + self.inflation_rate)
yearly_balances.append(balance)
# Early termination if portfolio is depleted
if balance <= 0:
yearly_balances += [0] * (self.years - year - 1)
break
simulations.append(yearly_balances)
return np.array(simulations)
def plot_simulation(self, simulations):
for sim in simulations:
plt.plot(sim, alpha=0.08, color='steelblue')
plt.title('Monte Carlo Portfolio Simulation')
plt.xlabel('Years')
plt.ylabel('Portfolio Value (Inflation Adjusted)')
plt.grid(True)
plt.show()
def summary_statistics(self, simulations):
ending_balances = simulations[:, -1]
print(f"\nSummary after {self.years} years across {self.num_simulations} simulations:")
print(f"Median Ending Balance: ${np.median(ending_balances):,.2f}")
print(f"10th Percentile: ${np.percentile(ending_balances, 10):,.2f}")
print(f"90th Percentile: ${np.percentile(ending_balances, 90):,.2f}")
print(f"Probability of Portfolio Depletion: {np.mean(ending_balances <= 0) * 100:.2f}%")