Part 3 of 7
🎯 Learning Objective
By the end of this section, you will have built a complete Monte Carlo risk simulation system in Python — from single-stock simulation to multi-asset portfolio VaR estimation.
| Purpose | Library |
|---|---|
| Numerical computation | numpy |
| Data manipulation | pandas |
| Visualization | matplotlib |
Stock prices are modeled using Geometric Brownian Motion (GBM):
Where: $S_t$ = stock price, $\mu$ = expected return (drift), $\sigma$ = volatility, $dW_t$ = Wiener process increment.
Where $Z \sim \mathcal{N}(0, 1)$ — a standard normal random variable.
💡 Intuition
Each future stock price is the current price, scaled by a deterministic growth factor ($\mu$) and a random shock ($\sigma \cdot Z$). The randomness creates the "fan" of possible future prices.
Pythonimport numpy as np
import matplotlib.pyplot as plt
# ─── Parameters ───
S0 = 100 # Initial stock price (₹)
mu = 0.08 # Expected annual return (8%)
sigma = 0.20 # Annual volatility (20%)
T = 1.0 # Time horizon (1 year)
N = 252 # Number of trading days
dt = T / N # Time step
n_simulations = 10000 # Number of Monte Carlo paths
# ─── Simulation ───
np.random.seed(42) # For reproducibility
# Generate random shocks
Z = np.random.standard_normal((n_simulations, N))
# Calculate daily returns
daily_returns = np.exp((mu - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * Z)
# Build price paths
price_paths = np.zeros((n_simulations, N + 1))
price_paths[:, 0] = S0
for t in range(1, N + 1):
price_paths[:, t] = price_paths[:, t-1] * daily_returns[:, t-1]
# ─── Visualization ───
plt.figure(figsize=(12, 6))
plt.plot(price_paths[:100].T, alpha=0.3, linewidth=0.8)
plt.title('Monte Carlo: 100 Stock Price Paths (GBM)', fontsize=14, fontweight='bold')
plt.xlabel('Trading Days')
plt.ylabel('Stock Price (₹)')
plt.axhline(y=S0, color='red', linestyle='--', label=f'Initial Price = ₹{S0}')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"After 1 year ({N} trading days):")
print(f" Mean final price: ₹{price_paths[:, -1].mean():.2f}")
print(f" Std of final price: ₹{price_paths[:, -1].std():.2f}")
print(f" Min final price: ₹{price_paths[:, -1].min():.2f}")
print(f" Max final price: ₹{price_paths[:, -1].max():.2f}")If σ = 0.40 instead of 0.20:
This is why volatility is the key input in risk models. Higher σ = wider distribution = more uncertainty = higher risk.
Python# ─── Calculate Returns ───
final_prices = price_paths[:, -1]
returns = (final_prices - S0) / S0 # Simple returns
log_returns = np.log(final_prices / S0) # Log returns
# ─── Visualization: Histogram of Returns ───
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# Simple Returns
axes[0].hist(returns, bins=80, density=True, alpha=0.7,
color='steelblue', edgecolor='white')
axes[0].axvline(x=returns.mean(), color='red', linestyle='--',
linewidth=2, label=f'Mean = {returns.mean():.2%}')
axes[0].axvline(x=0, color='black', linestyle='-', linewidth=1)
axes[0].set_title('Distribution of Annual Returns', fontweight='bold')
axes[0].set_xlabel('Return')
axes[0].set_ylabel('Density')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# Log Returns
axes[1].hist(log_returns, bins=80, density=True, alpha=0.7,
color='darkorange', edgecolor='white')
axes[1].axvline(x=log_returns.mean(), color='red', linestyle='--',
linewidth=2, label=f'Mean = {log_returns.mean():.2%}')
axes[1].axvline(x=0, color='black', linestyle='-', linewidth=1)
axes[1].set_title('Distribution of Log Returns', fontweight='bold')
axes[1].set_xlabel('Log Return')
axes[1].set_ylabel('Density')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()Three key reasons:
| Property | Simple Returns | Log Returns |
|---|---|---|
| Distribution | Right-skewed | Approximately normal |
| Additive over time? | ❌ No (must compound) | ✅ Yes |
| Can go below -100%? | ❌ No (bounded at -1) | ✅ Yes (−∞ to +∞) |
| Symmetry | Asymmetric | Symmetric |
Log returns are time-additive and approximately normally distributed — both critical properties for GBM modeling and Monte Carlo simulation.
Python# ─── VaR Calculation ───
portfolio_value = 1_000_000 # ₹10 Lakh portfolio
# Calculate P&L for each simulation
pnl = (final_prices / S0 - 1) * portfolio_value
# Compute VaR at different confidence levels
confidence_levels = [0.90, 0.95, 0.99]
print("=" * 60)
print(f"{'Monte Carlo VaR Analysis':^60}")
print(f"{'Portfolio Value: ₹10,00,000':^60}")
print(f"{'Simulations: ' + str(n_simulations):^60}")
print("=" * 60)
print(f"{'Confidence Level':<20} {'VaR (₹)':<15} {'CVaR (₹)':<15}")
print("-" * 60)
for cl in confidence_levels:
var = -np.percentile(pnl, (1 - cl) * 100)
cvar = -pnl[pnl <= -var].mean()
print(f"{cl:.0%}{'':<16} {var:>12,.0f} {cvar:>12,.0f}")
print("=" * 60)Python# ─── Visualization: VaR on Distribution ───
plt.figure(figsize=(12, 6))
plt.hist(pnl, bins=100, density=True, alpha=0.7,
color='steelblue', edgecolor='white')
# Mark VaR levels
var_95 = -np.percentile(pnl, 5)
var_99 = -np.percentile(pnl, 1)
cvar_95 = -pnl[pnl <= -var_95].mean()
plt.axvline(x=-var_95, color='orange', linestyle='--', linewidth=2,
label=f'95% VaR = ₹{var_95:,.0f}')
plt.axvline(x=-var_99, color='red', linestyle='--', linewidth=2,
label=f'99% VaR = ₹{var_99:,.0f}')
plt.axvline(x=-cvar_95, color='darkred', linestyle=':', linewidth=2,
label=f'95% CVaR = ₹{cvar_95:,.0f}')
plt.title('Portfolio P&L Distribution with VaR and CVaR',
fontsize=14, fontweight='bold')
plt.xlabel('Profit & Loss (₹)')
plt.ylabel('Density')
plt.legend(fontsize=10)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()Question: The 95% VaR is ₹2,12,000 on a ₹10 Lakh portfolio. What does this mean in practice?
Answer:
Risk manager's perspective: This portfolio needs at least ₹3,24,000 in risk capital reserves to cover 99% of annual scenarios.
For a portfolio of multiple assets, we account for correlations:
Where $\mathbf{w}$ = weight vector and $\boldsymbol{\Sigma}$ = covariance matrix.
Pythonimport pandas as pd
# ─── Define Multi-Asset Portfolio ───
assets = ['RELIANCE', 'TCS', 'HDFCBANK', 'INFY', 'ITC']
expected_returns = np.array([0.12, 0.15, 0.10, 0.14, 0.08])
volatilities = np.array([0.25, 0.22, 0.18, 0.24, 0.15])
correlation_matrix = np.array([
[1.00, 0.45, 0.50, 0.40, 0.30],
[0.45, 1.00, 0.35, 0.70, 0.20],
[0.50, 0.35, 1.00, 0.30, 0.40],
[0.40, 0.70, 0.30, 1.00, 0.25],
[0.30, 0.20, 0.40, 0.25, 1.00]
])
# Construct covariance matrix: Σ = D * C * D
D = np.diag(volatilities)
cov_matrix = D @ correlation_matrix @ D
# ─── Portfolio Weights ───
weights = np.array([0.25, 0.20, 0.25, 0.15, 0.15])
# ─── Portfolio Expected Return & Risk ───
portfolio_return = np.dot(weights, expected_returns)
portfolio_variance = weights @ cov_matrix @ weights
portfolio_volatility = np.sqrt(portfolio_variance)
print(f"Portfolio Expected Annual Return: {portfolio_return:.2%}")
print(f"Portfolio Annual Volatility: {portfolio_volatility:.2%}")
print(f"Sharpe Ratio (Rf=6%): {(portfolio_return - 0.06) / portfolio_volatility:.2f}")Python# ─── Correlated Multi-Asset Simulation ───
n_sims = 10000
n_days = 252
portfolio_value = 10_000_000 # ₹1 Crore
# Cholesky decomposition for correlated random variables
L = np.linalg.cholesky(cov_matrix)
daily_returns_all = np.zeros((n_sims, n_days, len(assets)))
for sim in range(n_sims):
Z = np.random.standard_normal((n_days, len(assets)))
correlated_Z = Z @ L.T
for i, (mu_i, sig_i) in enumerate(zip(expected_returns, volatilities)):
daily_returns_all[sim, :, i] = np.exp(
(mu_i - 0.5 * sig_i**2) * dt + correlated_Z[:, i] * np.sqrt(dt)
) - 1
# Calculate portfolio daily returns
portfolio_daily_returns = np.sum(daily_returns_all * weights, axis=2)
portfolio_values = portfolio_value * np.cumprod(1 + portfolio_daily_returns, axis=1)
final_portfolio_values = portfolio_values[:, -1]
portfolio_pnl = final_portfolio_values - portfolio_value
# ─── Portfolio VaR ───
print("\n" + "=" * 65)
print(f"{'Multi-Asset Portfolio Monte Carlo VaR':^65}")
print(f"{'Portfolio Value: ₹1,00,00,000 (₹1 Crore)':^65}")
print("=" * 65)
for cl in [0.90, 0.95, 0.99]:
var = -np.percentile(portfolio_pnl, (1 - cl) * 100)
cvar = -portfolio_pnl[portfolio_pnl <= -var].mean()
print(f" {cl:.0%} VaR: ₹{var:>12,.0f} | CVaR: ₹{cvar:>12,.0f}")
print("=" * 65)Problem: We need random numbers for 5 assets that are correlated with each other, not independent.
Solution — Cholesky Decomposition:
Why it works: $\text{Cov}(\tilde{Z}) = L \cdot \text{Cov}(Z) \cdot L^T = L \cdot I \cdot L^T = L L^T = \Sigma$ ✅
Without this, our simulation would treat all assets as independent — drastically underestimating (or overestimating) portfolio risk.
⚠️ Critical Teaching Point
Before proceeding to the quantum section, this is the most important concept to understand.
Python# ─── Demonstrate the Convergence Problem ───
import time
sim_counts = [100, 1000, 5000, 10000, 50000, 100000]
var_estimates = []
times = []
for n in sim_counts:
start = time.time()
Z = np.random.standard_normal(n)
sim_returns = np.exp((mu - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)
sim_pnl = (sim_returns - 1) * portfolio_value
var_est = -np.percentile(sim_pnl, 5)
elapsed = time.time() - start
var_estimates.append(var_est)
times.append(elapsed)
print("\nConvergence Analysis:")
print(f"{'Simulations':<15} {'95% VaR (₹)':<18} {'Time (ms)':<12} {'Relative Error'}")
print("-" * 60)
for n, var, t in zip(sim_counts, var_estimates, times):
rel_error = abs(var - var_estimates[-1]) / var_estimates[-1]
print(f"{n:<15,} {var:<18,.0f} {t*1000:<12.2f} {rel_error:<.4%}")💡 The Bridge to Quantum
"Classical Monte Carlo works. It's reliable. It's the industry standard. But look at the convergence table — to cut the error in half, you need 4× the simulations. To cut it by 10×, you need 100× the simulations."
| Metric | Classical MC | Quantum AE |
|---|---|---|
| Convergence rate | $O(1/\sqrt{N})$ | $O(1/N)$ |
| 10× better accuracy needs | 100× more sims | 10× more calls |
| For 100,000 positions overnight | Hours of compute | Minutes (future) |
"For a bank running 100,000 positions overnight, this $O(1/\sqrt{N})$ convergence creates a real computational wall. Now let's understand how quantum computing proposes to break through this wall."