Monte Carlo Simulation Explained

Monte Carlo simulation is one of the most powerful tools in quantitative finance. It works by running thousands of random scenarios to estimate the range of possible outcomes. If you cannot solve a problem analytically, simulate it.

Intuition: The Coin Flip Example

Suppose you flip a coin 100 times and bet $1 each time. You win $1 on heads, lose $1 on tails. What is the range of possible outcomes? We could do the maths, or we could just simulate it 10,000 times.

import numpy as np import matplotlib.pyplot as plt np.random.seed(42) num_simulations = 10000 num_flips = 100 # Simulate: +1 for heads, -1 for tails results = np.random.choice([1, -1], size=(num_simulations, num_flips)) # Cumulative sum for each simulation final_values = results.sum(axis=1) # Results print(f"Mean outcome: ${{final_values.mean():.2f}}") print(f"Std deviation: ${{final_values.std():.2f}}") print(f"Best case:  ${{final_values.max()}}") print(f"Worst case: ${{final_values.min()}}") print(f"Probability of profit: {{(final_values > 0).mean():.2%}}") # Histogram plt.figure(figsize=(10, 5)) plt.hist(final_values, bins=50, color="teal", edgecolor="black", alpha=0.7) plt.axvline(x=0, color="red", linewidth=2, linestyle="--") plt.title("Monte Carlo: 10,000 Simulations of 100 Coin Flips") plt.xlabel("Final P&L ($)") plt.ylabel("Frequency") plt.grid(True, alpha=0.3) plt.tight_layout() plt.show()

Financial Application: Portfolio Projection

The real power of Monte Carlo: projecting a portfolio over 30 years to answer retirement questions like “Will I have enough?”

import numpy as np import matplotlib.pyplot as plt # Parameters initial_investment = 100000   # $100,000 annual_return = 0.08          # 8% expected annual return annual_volatility = 0.16      # 16% annual volatility (typical for stocks) years = 30 trading_days = 252 num_simulations = 5000 # Daily parameters daily_return = annual_return / trading_days daily_vol = annual_volatility / np.sqrt(trading_days) # Simulate np.random.seed(42) total_days = years * trading_days random_returns = np.random.normal(daily_return, daily_vol, (num_simulations, total_days)) # Compound returns to get portfolio paths portfolio_paths = initial_investment * np.cumprod(1 + random_returns, axis=1) # Final values final_values = portfolio_paths[:, -1] # Key percentiles percentiles = [10, 25, 50, 75, 90] for p in percentiles: val = np.percentile(final_values, p) print(f"{{p}}th percentile: ${{val:,.0f}}") print(f"\nMean: ${{final_values.mean():,.0f}}") print(f"Probability of at least doubling: {{(final_values >= 2 * initial_investment).mean():.1%}}") # Plot a sample of paths plt.figure(figsize=(12, 6)) for i in range(100):  # Plot first 100 paths plt.plot(portfolio_paths[i], alpha=0.05, color="teal", linewidth=0.5) # Plot key percentiles time_axis = np.arange(total_days) for p in [10, 50, 90]: path_percentile = np.percentile(portfolio_paths, p, axis=0) plt.plot(time_axis, path_percentile, linewidth=2, label=f"{{p}}th percentile: ${{path_percentile[-1]:,.0f}}") plt.title(f"Monte Carlo: ${{initial_investment:,.0f}} Portfolio Over {{years}} Years") plt.xlabel("Trading Days") plt.ylabel("Portfolio Value ($)") plt.legend() plt.grid(True, alpha=0.3) plt.tight_layout() plt.show()

Adding Monthly Contributions

# More realistic: add $500/month to the portfolio monthly_contribution = 500 # Simulate again with contributions portfolio_paths_contrib = np.zeros((num_simulations, total_days)) portfolio_paths_contrib[:, 0] = initial_investment for day in range(1, total_days): # Add monthly contribution on the first trading day of each month (~every 21 days) contribution = monthly_contribution if day % 21 == 0 else 0 daily_r = np.random.normal(daily_return, daily_vol, num_simulations) portfolio_paths_contrib[:, day] = ( portfolio_paths_contrib[:, day - 1] * (1 + daily_r) + contribution ) final_with_contrib = portfolio_paths_contrib[:, -1] total_contributed = initial_investment + monthly_contribution * 12 * years print(f"Total money put in: ${{total_contributed:,.0f}}") print(f"Median final value: ${{np.percentile(final_with_contrib, 50):,.0f}}") print(f"90% confidence minimum: ${{np.percentile(final_with_contrib, 10):,.0f}}")