NumPy: Project — Monte Carlo Simulation
Last updated: 2026-08-26
1. Project: Monte Carlo Simulation
(1) Scenario
Alice needs to assess portfolio risk — Value at Risk (VaR) has no closed-form formula. Traditional variance-covariance methods assume normal distributions, severely underestimating tail risk. She uses Monte Carlo simulation to model 100,000 price paths.
(1) The Monte Carlo Method
Core idea: use many random samples to approximate real values. When problems are too complex for analytical solutions, simulation is the only way.
PYTHON
import numpy as np
rng = np.random.default_rng(42)
N = 100_000
x = rng.random(N)
y = rng.random(N)
inside = (x**2 + y**2) <= 1.0
pi_est = 4 * inside.mean()
print(f"Pi estimate: {pi_est:.6f}")
print(f"Actual pi: {np.pi:.6f}")
print(f"Error: {abs(pi_est - np.pi):.6f}")
(2) Estimating Pi
PYTHON
print("=== Monte Carlo Estimation of Pi ===\n")
for N in [1_000, 10_000, 100_000, 1_000_000]:
x = rng.random(N)
y = rng.random(N)
inside = (x**2 + y**2) <= 1.0
pi_est = 4 * inside.mean()
error = abs(pi_est - np.pi)
print(f"N={N:>9,d} Pi≈{pi_est:.6f} Error={error:.6f}")
(3) Random Walk
PYTHON
N = 1000
n_sim = 10_000
steps = rng.choice([-1, 1], size=(n_sim, N))
final_pos = steps.sum(axis=1)
print(f"Final position statistics ({n_sim} simulations):")
print(f" Mean: {final_pos.mean():.2f} (theory: 0)")
print(f" Std: {final_pos.std():.2f} (theory: {np.sqrt(N):.2f})")
print(f" Min: {final_pos.min()}")
print(f" Max: {final_pos.max()}")
(4) Option Pricing
PYTHON
S0 = 100
K = 105
r = 0.05
sigma = 0.20
T = 1.0
for n_sim in [10_000, 100_000, 1_000_000]:
z = rng.standard_normal(n_sim)
ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * z)
call_payoff = np.maximum(ST - K, 0)
put_payoff = np.maximum(K - ST, 0)
call_price = np.exp(-r * T) * call_payoff.mean()
put_price = np.exp(-r * T) * put_payoff.mean()
call_se = np.exp(-r * T) * call_payoff.std() / np.sqrt(n_sim)
print(f"N={n_sim:>9,d}")
print(f" Call: {call_price:.4f} ± {1.96*call_se:.4f}")
print(f" Put: {put_price:.4f}")
(5) Bootstrap Confidence Intervals
PYTHON
data = rng.exponential(scale=2.0, size=50)
n = len(data)
B = 10_000
idx = rng.integers(0, n, size=(B, n))
boot_samples = data[idx]
boot_means = boot_samples.mean(axis=1)
ci = np.percentile(boot_means, [2.5, 97.5])
print(f"Bootstrap 95% CI for mean: [{ci[0]:.3f}, {ci[1]:.3f}]")
print(f"Sample mean: {data.mean():.3f}")
▶ Example: Estimating Pi with Monte Carlo (Difficulty ⭐⭐)
PYTHON
import numpy as np
rng = np.random.default_rng(42)
N = 100_000
x = rng.random(N)
y = rng.random(N)
inside = (x**2 + y**2) <= 1.0
pi_est = 4 * inside.mean()
print(f"Estimated Pi: {pi_est:.6f}")
print(f"Actual Pi: {np.pi:.6f}")
print(f"Error: {abs(pi_est - np.pi):.6f}")
Output:
TEXT 📖 Display onlyEstimated Pi: 3.141920 Actual Pi: 3.141593 Error: 0.000327
▶ Example: Simulating a random walk (Difficulty ⭐⭐)
PYTHON
import numpy as np
rng = np.random.default_rng(42)
n_steps = 1000
n_simulations = 100
# All walks at once using vectorization
steps = rng.choice([-1, 1], size=(n_simulations, n_steps))
positions = np.cumsum(steps, axis=1)
final = positions[:, -1]
print(f"Mean final position: {final.mean():.2f}")
print(f"Std final position: {final.std():.2f}")
print(f"Min: {final.min()}, Max: {final.max()}")
Output:
TEXT 📖 Display onlyMean final position: -0.80 Std final position: 31.72 Min: -78, Max: 72
▶ Example: Bootstrap confidence intervals (Difficulty ⭐⭐⭐)
PYTHON
import numpy as np
rng = np.random.default_rng(42)
data = rng.exponential(scale=2.0, size=50)
n = len(data)
B = 10_000
indices = rng.integers(0, n, size=(B, n))
boot_means = data[indices].mean(axis=1)
ci = np.percentile(boot_means, [2.5, 97.5])
print(f"Sample mean: {data.mean():.3f}")
print(f"95% CI: [{ci[0]:.3f}, {ci[1]:.3f}]")
print(f"CI width: {ci[1] - ci[0]:.3f}")
Output:
TEXT 📖 Display onlySample mean: 2.165 95% CI: [1.738, 2.644] CI width: 0.906
❓ FAQ
Q How many Monte Carlo simulations do I need?
A Error ∝ 1/√N. 10,000 gives ~1% error, 1,000,000 gives ~0.1%. Financial pricing typically uses 1M, quick validation 10K. Run a small sample first to estimate the standard error, then scale up.
Q What's the difference between VaR and CVaR?
A VaR (Value at Risk) is the α-percentile loss threshold — "the worst 5% starts here." CVaR (Conditional VaR) is the average loss in the worst 5% — "when things go bad, how bad on average?" CVaR is always ≥ VaR.
Q Are Monte Carlo results reliable?
A Two caveats: (1) pseudo-random numbers are deterministic — different seeds give different results, always report confidence intervals; (2) model assumptions matter — if your GBM assumption is wrong, a million simulations won't save you.
📖 Summary
In this project you applied:
- Monte Carlo simulation: random sampling + aggregation
- Pi estimation: geometric probability with
rng.random - Random walks:
cumsumfor path simulation - Option pricing: vectorized payoff computation
- Bootstrap:
rng.choicewith replacement for confidence intervals - All powered by NumPy vectorization — no Python loops
📝 Exercises
-
Beginner (Difficulty ⭐): Modify the Pi estimation to use 1,000,000 points. How does the error change compared to 100,000?
-
Intermediate (Difficulty ⭐⭐): Implement a 2D random walk. Each step moves in one of 4 directions (up/down/left/right). Simulate 10,000 paths of 1000 steps and compute the mean final distance from the origin.
-
Advanced (Difficulty ⭐⭐⭐): Price an Asian option (payoff depends on average price over the path) using Monte Carlo. Compare the price to a European option with the same parameters. Explain why the Asian option is cheaper.
